An accurate low-redshift measurement of the cosmic neutral hydrogen density

TL;DR

This study uses spectral stacking of radio data from SDSS galaxies to precisely measure the cosmic neutral hydrogen density at low redshift, confirming its stability over recent cosmic times.

Contribution

It provides a new, robust estimate of the cosmic HI mass density at low redshift using a carefully calibrated stacking technique and addresses systematic uncertainties.

Findings

Measured $ ho_{HI}$ at $z hickapprox 0.07$ with high precision.

Found no significant evolution of HI content at low redshift.

Results are consistent with previous surveys and stacking experiments.

Abstract

Using a spectral stacking technique, we measure the neutral hydrogen (HI) properties of a sample of galaxies at $z < 0.11$ across 35 pointings of the Westerbork Synthesis Radio Telescope (WSRT). The radio data contains 1,895 galaxies with redshifts and positions known from the Sloan Digital Sky Survey (SDSS). We carefully quantified the effects of sample bias, aperture used to extract spectra, sidelobes and weighting technique and use our data to provide a new estimate for the cosmic HI mass density. We find a cosmic HI mass density of $\Omega_{\rm HI} = (4.02 \pm 0.26)\times 10^{-4} h_{70}^{-1}$ at $\langle z\rangle = 0.066$, consistent with measurements from blind HI surveys and other HI stacking experiments at low redshifts. The combination of the small interferometer beam size and the large survey volume makes our result highly robust against systematic effects due to confusion at small scales and cosmic variance at large scales. Splitting into three sub-samples with $\langle z\rangle$ = 0.038, 0.067 and 0.093 shows no significant evolution of the HI gas content at low redshift.

Tables

Table 1. Table 1: Mean HI mass M HI subscript 𝑀 HI M_{\rm HI} and mean r 𝑟 r -band HI mass-to-light ratio M HI / L r subscript 𝑀 HI subscript 𝐿 𝑟 M_{\rm HI}/L_{r} for the magnitude-limited sample using different values of the weight parameter γ 𝛾 \gamma . The weighted mean redshift is also given.

$\gamma$	<z>	<$M_{\mathrm{HI}}$>	$S/N\left(M_{\mathrm{HI}}\right)$	<$M_{\mathrm{HI}}/L_r$>	$S/N\left(M_{\mathrm{HI}}/L\right)$
		(${10}^9{h}_{70}^{-2}{M}_{\odot }$)		($M_{\odot }/L_{\odot }$)
0	0.062	2.75 $\pm$ 0.20	13.8	0.28 $\pm$ 0.03	10.4
1	0.051	2.34 $\pm$ 0.14	16.3	0.31 $\pm$ 0.02	15.5
2	0.041	1.90 $\pm$ 0.13	14.4	0.34 $\pm$ 0.04	7.7
3	0.032	1.55 $\pm$ 0.20	7.6	0.35 $\pm$ 0.08	4.6
4	0.025	1.33 $\pm$ 0.27	4.8	0.33 $\pm$ 0.07	4.5

Table 2. Table 2: Average HI mass M HI subscript 𝑀 HI M_{\rm HI} , mass-to-light ratio M HI / L subscript 𝑀 HI 𝐿 M_{\rm HI}/L , and corresponding signal-to-noise ratio S / N 𝑆 𝑁 S/N as a function of aperture size.

Aperture	<$M_{\mathrm{HI}}$>	$S/N\left(M_{\mathrm{HI}}\right)$	<$M_{\mathrm{HI}}/L$>	$S/N\left(M_{\mathrm{HI}}/L\right)$
(kpc)	(${10}^9{h}_{70}^{-2}{M}_{\odot }$)		($M_{\odot }{L}_{\odot }^{-1}$)
10	1.31 $\pm$ 0.08	15.5	0.19 $\pm$ 0.03	7.5
15	1.54 $\pm$ 0.12	12.9	0.22 $\pm$ 0.02	9.1
20	1.78 $\pm$ 0.12	14.9	0.25 $\pm$ 0.03	9.4
25	2.00 $\pm$ 0.15	13.4	0.28 $\pm$ 0.02	12.1
30	2.18 $\pm$ 0.14	15.1	0.30 $\pm$ 0.03	11.0
35	2.34 $\pm$ 0.14	16.3	0.31 $\pm$ 0.02	15.5
40	2.46 $\pm$ 0.16	15.0	0.32 $\pm$ 0.02	13.2
45	2.59 $\pm$ 0.13	19.6	0.33 $\pm$ 0.02	14.6
50	2.70 $\pm$ 0.18	15.0	0.35 $\pm$ 0.04	8.1
55	2.78 $\pm$ 0.21	13.3	0.36 $\pm$ 0.06	6.4
60	2.86 $\pm$ 0.25	11.6	0.38 $\pm$ 0.04	10.5
65	2.91 $\pm$ 0.19	15.6	0.39 $\pm$ 0.04	11.1
70	2.97 $\pm$ 0.20	15.2	0.40 $\pm$ 0.04	10.5
75	3.03 $\pm$ 0.24	12.7	0.41 $\pm$ 0.06	6.4
80	3.07 $\pm$ 0.17	18.5	0.40 $\pm$ 0.06	6.7
85	3.20 $\pm$ 0.23	13.6	0.44 $\pm$ 0.08	5.3

Table 3. Table 3: Measurements of ⟨ M HI ⟩ delimited-⟨⟩ subscript 𝑀 HI \langle M_{\rm HI}\rangle and ⟨ M HI / L ⟩ delimited-⟨⟩ subscript 𝑀 HI 𝐿 \langle M_{\rm HI}/L\rangle from stacking galaxies in the mock catalog. The superscript ‘m’ and ‘v’ refer to the magnitude-limited sample and volume-limited sub-sample respectively.

	No confusion	Confused with	Confused with
		sample galaxies	all galaxies
$⟨M_{\mathrm{HI}}^m⟩({10}^9{h}_{70}^{-2}$ M⊙)	2.757	2.794	2.816
$⟨M_{\mathrm{HI}}^m/L⟩$(M${}_{\odot }{}^{}/$L⊙)	0.289	0.293	0.294
$⟨M_{\mathrm{HI}}^v⟩({10}^9{h}_{70}^{-2}$ M⊙)	1.674	1.706	1.708
$⟨M_{\mathrm{HI}}^v/L⟩$(M${}_{\odot }{}^{}/$L⊙)	0.560	0.570	0.572

Table 4. Table 4: The results of stacking with the original catalogue, the simulated sky and the convolved sky. For the latter two, we use a aperture with a radius of 35 kpc to extract the spectra.

Data source	Aperture	Stacked integral
	(kpc)	(${10}^9{h}_{70}^{-2}$ M⊙)
$S^3$-SAX catalogue	–	3.013
Confused sky	35	3.021
Convolved sky	35	2.962

Table 5. Table 5: Results from stacking HI spectra in the 35 individual pointings observed with the WRST. Because of the smaller sample, the effective volume reduces and cosmic variance increases, the errors are larger compared with the results from the stacking with total sample(section  3.2 ). We also show the statistical errors for stacked mass and mass-to-light spectra as noise stat,m and noise stat,m/l .

Pointing	Position	$N$	<$z$>	<$M_{\mathrm{HI}}$>	noisestat,m	<$M_{\mathrm{HI}}/L$>	noisestat,m/l	Obs time
	(J2000)			(${10}^9{h}_{70}^{-2}$ M⊙)	(${10}^9{h}_{70}^{-2}$ M⊙)	(M${}_{\odot }{}^{}/$L⊙)	(M${}_{\odot }{}^{}/$L⊙)	(hrs)
1	22:27:00 +13:37:48	36	0.092	5.15 $\pm$ 3.96	0.41	1.33 $\pm$ 1.11	0.11	12.0
2	22:37:48 +14:18:36	66	0.074	3.43 $\pm$ 0.57	0.24	0.22 $\pm$ 0.08	0.03	12.0
3	22:57:50 +13:03:36	45	0.057	1.80 $\pm$ 0.53	0.11	0.70 $\pm$ 0.31	0.04	12.0
4	23:12:58 +13:56:24	71	0.066	5.38 $\pm$ 0.84	0.20	0.75 $\pm$ 0.21	0.04	11.5
5	23:14:24 +14:39:00	49	0.074	4.94 $\pm$ 1.42	0.30	0.45 $\pm$ 0.10	0.03	11.0
6	23:24:54 +15:18:00	70	0.056	3.18 $\pm$ 0.55	0.11	0.35 $\pm$ 0.10	0.02	10.7
7	23:43:23 +14:16:08	36	0.073	11.13 $\pm$ 4.46	0.59	0.76 $\pm$ 0.37	0.05	9.8
8	23:51:36 +14:06:00	46	0.078	5.17 $\pm$ 1.15	0.27	0.53 $\pm$ 0.10	0.04	8.8
9	02:03:18 +13:51:00	31	0.063	2.28 $\pm$ 1.18	0.31	0.18 $\pm$ 0.08	0.03	11.3
10	22:12:29 +12:20:24	31	0.067	1.44 $\pm$ 0.48	0.17	0.44 $\pm$ 0.17	0.04	12.0
11	22:14:38 +13:52:12	81	0.044	0.33 $\pm$ 0.14	0.05	–	–	9.7
12	22:33:18 +13:11:02	35	0.089	3.98 $\pm$ 1.41	0.62	0.32 $\pm$ 0.09	0.04	12.0
13	22:39:00 +13:26:24	57	0.079	1.58 $\pm$ 1.07	0.28	0.60 $\pm$ 0.74	0.06	12.0
14	23:18:18 +14:55:12	39	0.081	3.41 $\pm$ 1.22	0.44	0.30 $\pm$ 0.10	0.08	10.7
15	23:26:24 +14:03:00	55	0.054	2.11 $\pm$ 0.49	0.16	0.43 $\pm$ 0.12	0.05	10.3
16	23:38:06 +15:45:43	60	0.066	2.53 $\pm$ 0.82	0.21	0.19 $\pm$ 0.08	0.03	8.6
17	23:45:36 +15:22:12	26	0.087	–	–	0.06 $\pm$ 0.21	0.08	9.3
18	23:56:53 +13:57:00	27	0.067	8.92 $\pm$ 1.93	0.29	0.70 $\pm$ 0.14	0.04	12.0
19	00:00:36 +15:24:36	28	0.077	1.60 $\pm$ 3.55	0.38	0.05 $\pm$ 0.16	0.02	5.4
20	00:06:00 +15:43:48	36	0.069	5.91 $\pm$ 2.71	0.29	0.32 $\pm$ 0.19	0.03	10.0
21	00:24:00 +14:12:00	18	0.060	2.03 $\pm$ 1.79	0.21	0.15 $\pm$ 0.34	0.06	6.1
22	00:43:01 +15:18:00	74	0.080	1.34 $\pm$ 0.71	0.17	0.10 $\pm$ 0.04	0.01	10.8
23	01:10:03 +13:59:49	91	0.061	0.22 $\pm$ 0.68	0.09	0.07 $\pm$ 0.09	0.01	12.0
24	01:11:28 +15:06:00	63	0.055	3.96 $\pm$ 1.28	0.20	0.47 $\pm$ 0.14	0.03	12.0
25	01:15:00 +14:28:48	66	0.064	3.94 $\pm$ 0.60	0.23	0.46 $\pm$ 0.11	0.03	4.6
26	01:55:48 +14:45:07	75	0.068	3.10 $\pm$ 0.69	0.23	0.64 $\pm$ 0.17	0.03	10.1
27	01:57:11 +13:09:00	51	0.057	3.93 $\pm$ 0.97	0.20	0.43 $\pm$ 0.09	0.04	5.7
28	02:12:00 +14:02:24	40	0.048	2.53 $\pm$ 0.84	0.11	0.84 $\pm$ 0.43	0.04	10.1
29	00:00:36 +14:33:00	64	0.086	1.89 $\pm$ 1.43	0.43	0.18 $\pm$ 0.08	0.03	9.4
30	00:30:36 +14:52:12	33	0.074	2.29 $\pm$ 1.68	0.37	0.09 $\pm$ 0.13	0.04	7.8
31	00:58:01 +14:50:24	54	0.074	4.21 $\pm$ 1.19	0.35	0.26 $\pm$ 0.08	0.03	9.5
32	01:19:48 +14:45:40	57	0.050	1.43 $\pm$ 0.36	0.14	0.36 $\pm$ 0.24	0.03	9.8
33	01:46:30 +13:51:00	42	0.062	2.88 $\pm$ 1.04	0.25	0.35 $\pm$ 0.13	0.04	12.0
34	01:49:26 +13:51:00	88	0.062	3.37 $\pm$ 0.52	0.14	0.41 $\pm$ 0.08	0.03	10.8
35	23:24:18 +14:40:48	154	0.052	0.51 $\pm$ 0.30	0.11	0.08 $\pm$ 0.05	0.03	9.4

Table 6. Table 6: The results of stacking the sample in different redshift bins.

$⟨z⟩$	redshift range	$N_g$	<$M_{\mathrm{HI}}$>	<$M_{\mathrm{HI}}/L$>
			(${10}^9{h}_{70}^{-2}$M⊙)	(M${}_{\odot }{}^{}/$L⊙)
0.024 $\pm$ 0.004	0.01 - 0.03	155	0.90 $\pm$ 0.12	0.38 $\pm$ 0.11
0.041 $\pm$ 0.004	0.03 - 0.05	439	2.34 $\pm$ 0.22	0.37 $\pm$ 0.05
0.062 $\pm$ 0.005	0.05 - 0.07	453	2.87 $\pm$ 0.31	0.25 $\pm$ 0.03
0.080 $\pm$ 0.005	0.07 - 0.09	448	2.97 $\pm$ 0.55	0.21 $\pm$ 0.04
0.097 $\pm$ 0.005	0.09 - 0.11	400	4.45 $\pm$ 0.94	0.22 $\pm$ 0.05

Table 7. Table 7: The measurement of Ω HI subscript Ω HI \Omega_{\rm HI} with different methods. The K c subscript 𝐾 𝑐 K_{c} refers to the correction factor for confusion.

Method	Formula	${\Omega }_{\mathrm{HI}}\left({10}^{-4}{h}_{70}^{-1}\right)$
Measured	$⟨\frac{M_{\mathrm{HI}}}{L}⟩\times {\rho }_L$	2.96 $\pm$ 0.19
Luminosity bias corrected	$C1\times ⟨\frac{M_{\mathrm{HI}}}{L}⟩\times {\rho }_L$	4.09 $\pm$ 0.27
Confusion corrected	$K_c^{-1}\times C1\times ⟨\frac{M_{\mathrm{HI}}}{L}⟩\times {\rho }_L$	4.02 $\pm$ 0.26

Table 8. Table 8: The cosmic HI density Ω HI subscript Ω HI \Omega_{\rm HI} in different redshift bins. The confusion correction has been applied.

$⟨z⟩$	$N_g$	<$M_{\mathrm{HI}}$>	<$M_{\mathrm{HI}}/L$>	C1	${\Omega }_{\mathrm{HI}}$
		(${10}^9{h}_{70}^{-2}$M⊙)	(M${}_{\odot }{}^{}/$L⊙)		(${10}^{-4}{h}_{70}^{-1}$)
0.038 $\pm$ 0.009	634	1.73 $\pm$ 0.17	0.38 $\pm$ 0.06	1.23	3.92 $\pm$ 0.63
0.067 $\pm$ 0.007	637	2.83 $\pm$ 0.34	0.21 $\pm$ 0.03	2.13	3.97 $\pm$ 0.61
0.093 $\pm$ 0.007	621	3.94 $\pm$ 0.63	0.22 $\pm$ 0.02	1.92	3.99 $\pm$ 0.36

Equations

$$S_{\nu }=\frac{{\Sigma }_x{\Sigma }_y{S}_{\nu }(x,y)}{{\Sigma }_x{\Sigma }_{y}B(x,y)},$$

$$S_{{\nu }_{\mathrm{res}}}=\frac{S_{{\nu }_{\mathrm{obs}}}}{(1+z)}.$$

$$m_{\mathrm{HI}}(\nu )=4.98\times 10^7{S}_{\nu }{D}_L^2{f}^{-1},$$

$$w_i=f^2{D}_L^{-\gamma }{\sigma }^{-2},$$

$$⟨m_{\mathrm{HI}}(\nu )⟩=\frac{{\sum }_{i=1}^n{w}_i{m}_{\mathrm{HI},\mathrm{i}}}{{\sum }_{i=1}^n{w}_i}.$$

$$M_{\mathrm{HI}}={\int }_{{\nu }_0-\Delta \nu }^{{\nu }_0+\Delta \nu }⟨m_{\mathrm{HI}}(\nu )⟩d\nu ,$$

$${\sigma }^2(⟨m_{\mathrm{HI}}⟩)=\frac{19}{20}\sum _{j=1}^{20}(⟨m_{\mathrm{HI}}⟩-⟨m_{\mathrm{HI}}^j⟩{)}^2,$$

$$S_i=S_s+{\Sigma }_c{f}_{1;c}{f}_{2:c}{S}_c,$$

$$C1=\frac{⟨M_{\mathrm{HI}}/L{⟩}_{\mathrm{all}}}{⟨M_{\mathrm{HI}}/L{⟩}_{\mathrm{obs}}}=\frac{{\sum }_{i=1}^N{w}_i}{{\sum }_{i=1}^N{w}_i(L_i/L_{\ast }{)}^{\beta }}\times \frac{\int (L/L_{\ast }{)}^{\beta }L{\varphi }_L(L)dL}{\int L{\varphi }_L(L)dL},$$

$${\varphi }_L(L)dL={\varphi }^{\ast }{\left(\frac{L}{L_{\ast }}\right)}^{\alpha }\exp \left(-\frac{L}{L_{\ast }}\right)\frac{dL}{L_{\ast }},$$

$${\rho }_{\mathrm{HI}}=C1\times ⟨\frac{M_{\mathrm{HI}}}{L}⟩\times {\rho }_L.$$

$${\Omega }_{\mathrm{HI}}=\frac{{\rho }_{\mathrm{HI}}}{{\rho }_{\mathrm{c},0}}=(4.02\pm 0.26)\times 10^{-4}{h}_{70}^{-1}.$$

$${\Omega }_{\mathrm{HI}}^{\mathrm{v}}=(3.50\pm 0.90)\times 10^{-4}{h}_{70}^{-1},0$$

$${\rho }_L(z)={\rho }_L(z_0)10^{0.4(P+Q)(z-z_0)},$$

$$\alpha (z)=\alpha (z_0),$$

$$M^{\star }(z)=M^{\star }(z_0)-Q(z-z_0),$$

$${\phi }^{\star }(z)={\phi }^{\star }(0)10^{0.4Pz},$$

$${\rho }_{\mathrm{HI}}$$

$$\approx$$

$${\Omega }_{\mathrm{HI}}=(4.01\pm 0.30)\times 10^{-4}{h}_{70}^{-1}.$$

Full text

††pubyear: 2018††pagerange: An accurate low-redshift measurement of the cosmic neutral hydrogen density–A

1 Introduction

To fully understand the formation and evolution of galaxies, it is important to study the accretion of gas from the intergalactic medium (IGM), galaxy mergers and galaxy interaction, and the depletion of gas through galactic fountains and outflow processes (Kereš et al., 2005; Sancisi et al., 2008; Marinacci et al., 2010). Cool gas drives star formation in galaxies as shown by the correlation between star-formation surface density ($\Sigma_{\rm SFR}$) and HI surface density ($\Sigma_{\rm HI}$) (Schmidt, 1959; Kennicutt, 1998), and the even tighter correlation with molecular hydrogen surface density ($\Sigma_{\rm H_{2}}$) (Bigiel et al., 2008; Schruba et al., 2010). Whilst the latter provides evidence for the important role of molecular clouds in controlling star formation (Solomon & Vanden Bout, 2005), because of the relatively short gas consumption time scales, it is the large-scale net inflow and condensation of cool gas which eventually forms the massive molecular clouds prior to star formation. Therefore study of both the atomic and molecular phases of cool gas in galaxies is crucial for the understanding of their star formation history.

There are a number of observation techniques we can use to measure HI gas content. At high redshifts, the damped Lyman-$\alpha$ (DLA) systems seem to indicate large reservoirs of HI whose column density can be deduced from DLA absorption profiles, thereby allowing determination of the cosmic HI mass density. At $z>1.65$, many DLA surveys have therefore been used to measure the cosmic HI gas density (Lanzetta et al., 1991; Prochaska et al., 2005; Noterdaeme et al., 2009; Noterdaeme et al., 2012; Songaila & Cowie, 2010; Zafar et al., 2013; Crighton et al., 2015; Neeleman et al., 2016; Bird et al., 2017). Their results show a significant evolution of HI gas content over cosmic time and that there is more HI gas at high redshifts. At $z<1.65$, Lyman-$\alpha$ absorption is only detected at ultraviolet (UV) wavelengths, so can only be observed with space-based telescopes. Rao et al. (2006); Rao et al. (2017) have identified candidate DLA systems through their metal absorption lines in the redshift range $0.11<z<1.65$. Their results indicate no clear evolution of cosmic HI gas density. However, the low incidence of DLAs per unit redshift at intermediate redshifts give rise to significant statistical uncertainties.

In the local Universe, the HI content is conveniently measured through the direct detection of the 21-cm hyperfine emission line. The large instantaneous field of view provided by modern multibeam receivers has made blind, large-area HI surveys possible. The HI Parkes All-Sky Survey (HIPASS) (Barnes et al., 2001) has detected HI emission from 5317 galaxies at $0<z<0.04$ over a sky area of 21,341 deg2 (Meyer et al., 2004; Wong et al., 2006), and the Arecibo Legacy Fast ALFA (ALFALFA) survey (Giovanelli et al., 2005) has detected $\sim$ 31500 galaxies out to $z=0.06$ over a sky area of approximately 7000 deg2 (Haynes et al., 2018). These large-area surveys allow for accurate measurement of the local HI mass function and the cosmic HI gas density. The measurements of HI density from these surveys are reasonably consistent with each other (Zwaan et al., 2005; Martin et al., 2010; Jones et al., 2018). However, directly measuring 21-cm emission of more distant individual galaxies is difficult with the current generation of single-dish radio telescopes, so this approach is limited to low redshift.

Individual deep 21-cm pointings have proven the feasibility of detecting HI galaxies outside the local Universe and up to $z\approx 0.3$ (Catinella et al., 2008; Zwaan et al., 2001; Verheijen et al., 2007; Fernández et al., 2016). However, in order to increase the chance of detection, the observed areas are often pre-selected. For example, Catinella & Cortese (2015) detected 39 galaxies up to $z=0.25$ with the 305-m Arecibo telescope, selecting them by presence of H$\alpha$ emission, disk morphology and isolation. Zwaan et al. (2001) and Verheijen et al. (2007) targeted galaxies in clusters at $z\approx 0.2$ with the Westerbork Synthesis Radio Telescope (WSRT). These samples are biased towards bright galaxies with high optical surface brightness, or in dense regions.

However, blind surveys to higher redshifts are time consuming. For example, the Arecibo Ultra Deep Survey (AUDS) (Freudling et al., 2011; Hoppmann et al., 2015) has so-far detected 103 galaxies with 400 hrs of integration time in the redshift range of $0<z<0.16$. The Cosmological Evolution Survey (COSMOS) HI Large Extragalactic Survey (CHILES) over the redshift range $z=0$ - 0.45 (Fernández et al., 2013, 2016) will be able to detect up to 300 galaxies with 1000 hours of observation time on the Very Large Array (VLA). However, even with such large integration times, these surveys have been limited to very small sky areas (1.35 deg2 for AUDS and 0.3 deg2 for CHILES), resulting in small effective volumes and large cosmic variance.

Next generation telescopes SKA pathfinder such as Australian Square Kilometre Array Pathfinder (ASKAP) (Johnston et al., 2008; Meyer, 2009), MeerKAT (Holwerda et al., 2012), Five-hundred-meter Aperture Spherical radio Telescope (FAST) (Nan et al., 2011; Duffy et al., 2008) and WSRT/Aperture Tile in Focus (APERTIF) (Oosterloo et al., 2009) will enable large-area surveys to significant depths. But less direct methods for measuring HI gas content at higher redshifts are also available using the technique of spectral stacking (Chengalur et al., 2001). The technique combines a large number of rest-frame spectra extracted from the radio data with redshifts and positions from optical catalogues. This allows the noise to be averaged down, and recovers a more significant spectral line signal, but averaged over a large sample of galaxies. By potentially accessing a larger number of galaxies, HI stacking can provide significantly large volumes, and much smaller cosmic variance.

Studies using the spectral stacking technique for galaxies outside the local Universe include those of Verheijen et al. (2007) and Lah et al. (2009) who examined galaxies in cluster environments out to $z=0.37$. Other observations have been used to study the properties of nearby galaxies, for example the relation between the HI content of a galaxy and its bulge (Fabello et al., 2011b) and correlations between the HI content, stellar mass and environment (Fabello et al., 2012; Brown et al., 2015, 2018), as well as the influence of an active galactic nucleus (AGN) (Fabello et al., 2011a; Geréb et al., 2013). The first attempt to use stacking to calculate the cosmic HI gas density $\rm\Omega_{HI}$, was presented by Lah et al. (2007) in the redshift range $0.218<z<0.253$ using the Giant Metrewave Radio Telescope (GMRT). A more recent HI stacking experiment was carried out by Delhaize et al. (2013) using HIPASS data and new observations from the Parkes telescope combined them with $\sim 18,300$ redshifts from the Two-Degree Field Galaxy Redshift Survey (2dFGRS) to obtain high signal-to-noise ratio detections out to a redshift of $z=0.13$.

Rhee et al. (2013) used data from WSRT and stacked a significantly smaller sample of 59 galaxies at $z\approx 0.1$ and 96 galaxies at $z\approx 0.2$. Rhee et al. (2016) cross-matched the zCOSMOS-bright catalogue with data from GMRT, obtaining a 474 galaxy sample at $z\approx 0.37$. With the stacking technique, they made a 3$\sigma$ detection of average HI mass. Rhee et al. (2018) used observations made with the GMRT to probe the HI gas content of 165 field galaxies in the VIMOS VLT Deep Survey (VVDS) 14h field at $z\approx 0.32$, resulting in a measurement of HI mass with a significance of 2.8$\sigma$. Kanekar et al. (2016) used the GMRT to stack HI emission from massive star-forming galaxies at $z\approx 1.18\textendash-1.34$, the highest redshift at which stacking has been attempted.

The technique of ‘intensity mapping’ can also be used to extend the HI survey limit to higher redshifts. Similar to stacking, this involves measuring the cross-power between radio and optical surveys (Pen et al., 2009), but uses the bulk emission fluctuations due to galaxy clustering over the surveyed region instead of individual galaxies. Observations conducted with the Green Bank Telescope (Chang et al., 2010; Masui et al., 2013), spanning the redshift range $0.6<z<1$ have highlighted the potential power of this technique. However, the accuracy of cosmic HI density measurements remains low, and there is a dependence on simulations of the wavelength-dependent bias of galaxies at optical and radio wavelengths (Wolz et al., 2017).

In this paper, we foreshadow some of the techniques which will be utilised in the future SKA pathfinder surveys to bridge the redshift gap $0.2<z<1.65$. We achieve this by using an interferometer in order to reduce problems arising from confusion that affect single-dish data. But we also cover a wide field of view by using multiple pointing centres in order to reduce cosmic variance, which has otherwise affected deep interferometer surveys. We obtain the radio data from WSRT (Geréb et al., 2015) and use a corresponding optical catalog from SDSS (York et al., 2000) containing 1895 galaxies within the sampled redshift range. Sample selection is not biased by environment, star formation, or any particular physical characteristic other than the optical magnitude limits of the SDSS.

Section 2 presents the observational data used in this paper. In Section 3 we present the spectral extraction and stacking methodology. In Section 4 we measure average HI mass and HI mass-to-light ratio for the sample, and various sub-samples in redshift and luminosity. In Section 5 we describe our measurement of $\Omega_{HI}$ and compare with existing results in the literature. Throughout this paper we use H${}_{\circ}=70$ km s*-1* Mpc*-1*, $\Omega_{\rm m}=0.3$ and $\Omega_{\Lambda}=0.7$.

2 Data

The HI observations were made using the Westerbork Synthesis Radio Telescope (WSRT). Thirty six pointing positions were selected according to the overall WSRT schedule with the only main constraint being that the sky overlap with footprint of the Galaxy Evolution eXplorer (GALEX) survey (Martin et al., 2005) and Sloan Digital Sky Survey (SDSS) South Galactic Cap region ($\rm 21^{h}<RA<2^{h}$ and $\rm 10^{\circ}<Dec<16^{\circ}$). 351 hours of observation time were used to observe the region, with each pointing observed for between 5 hr and 12 hr. Data from one of the pointings were discarded due to bad data quality. The sky region covered by the remaining 35 pointings is shown in Figure 1.

The half-power beam width (HPBW) of the WSRT is 35 arcmin at the observing frequency, and the average synthesized beam size is $108\arcsec\times 22\arcsec$. Figure 2 shows a histogram of major axis, minor axis and position angles of the synthesized beams for the 35 pointings. The data were reduced and self-calibrated using the radio astronomy data reduction package miriad (Sault et al., 1995). The data were flagged to reduce the contamination by radio frequency interference (RFI).

The reduced data cubes have a size of 601$\times$601 pixels with a pixel size of $3\arcsec\times 3\arcsec$. The data consist of 8$\times$20 MHz bands, each with 128 channels and two polarisations. Each channel is 0.15625 MHz wide, corresponding to $\sim 33$ km s*-1* at $z=0$ and $\sim 37$ km s*-1* at $z=0.11$. The rms was typically 0.2 mJy beam*-1* per 0.15625 MHz channel for each field, independent of frequency. Each frequency band overlaps by 3 MHz resulting in an overall frequency range of 1.406 GHz to 1.268 GHz, corresponding to a redshift range of $0.01<z<0.12$. However, due to stronger RFI at higher redshift we set an upper redshift limit of $z=0.11$.

Accurate measurements of redshift and spacial positions are indispensable for stacking. We use SDSS DR9 as the optical catalogue for our stacking analysis. SDSS has a typical redshift error of $\sim 60$ km s*-1* and a spectral density of 60 - 100 deg*-2* ($z<0.12$) in the region we selected for the HI observations. With the target selection algorithm described in Strauss et al. (2002), the SDSS sample has a completeness which exceeds 99$\%$ (excluding fibre collisions). The sample appears to be complete for a star formation rate above 10*-2M☉* yr*-1* for $z<0.06$. The luminosities used in this paper are calculated from the SDSS $r$-band magnitudes, applying $k$-corrections (Chilingarian & Zolotukhin, 2012).

By cross-matching our radio data with the SDSS catalog, we obtain a sample of 1895 galaxies spanning the redshift range $0.01<z<0.11$ (Figure 3) and within the radius of the pointings at which the normalized primary beam response drops to 0.1. We refer to this as the magnitude-limited sample, only including galaxies with $r$-band magnitude brighter than 17.77. It has a mean redshift of $\langle z\rangle=0.066$. To measure the HI density with a sample less biased by magnitude, we also created a volume-limited sample with $z\leq 0.0285$, which has 149 galaxies in total and a mean redshift of $\langle z\rangle$ = 0.024. The volume-limited sample is complete for $r$-band luminosities $>10^{8.68}$L*⊙*. Figure 4 shows the $r$-band luminosity distribution as a function of redshift with the volume-limited sub-sample highlighted.

3 Stacking Analysis

3.1 HI Mass Spectra

The stacking technique used in this paper is similar to that described in Geréb et al. (2013). Spectra were extracted from the data cubes over an extended region around the SDSS position. After extensive tests, we find the region with aperture radius of 35kpc gives best stacking results(see Section 3.3). The spatially-integrated spectrum was calculated from:

[TABLE]

where $S_{\nu}(x,y)$ is the flux density at pixel position $(x,y)$ and $B(x,y)$ is the normalized synthesized beam response (centred on the SDSS position) at the same pixel position. After this, a second-order baseline was fitted to remove residual continuum (excluding a velocity range of 500 km s*-1* around the expected spectral location of the SDSS galaxy), and the spectra were de-redshifted. The barycentric frequency is converted from the observed to the rest frame by $\nu_{\rm res}=\nu_{\rm obs}(1+z)$. As the channel width is also broadened in this process, HI flux density is conserved by applying the corresponding correction:

[TABLE]

After shifting to the rest frame, the flux spectra were converted into mass spectra using the following relation:

[TABLE]

where $S_{\nu}$ is the de-redshifted HI flux density in Jy, $D_{L}$ is the luminosity distance in Mpc, $f$ is the normalised primary beam response, and $m_{\rm HI}$ is in units of M*☉* MHz*-1*.

We introduce a weight factor which depends on the primary beam response $f$, the luminosity distance $D_{L}$, as well as the rms noise of the flux density spectra $\sigma$. The weight of $i$-th galaxy is expressed as:

[TABLE]

where large values of $\gamma$ give larger weight to nearby galaxies, and small values give more weight to distant galaxies. The effect of the weight factor on the results is considered later. The averaged final stacked spectrum is obtained from:

[TABLE]

The integrated HI mass of a stack, or $\langle M_{\rm HI}\rangle$, is then defined as the integral along the frequency axis over the mass spectrum:

[TABLE]

where $\nu_{0}$ refers to 1420.406 MHz and $\Delta\nu$ is large enough to capture all flux from the stack (we will later use $\Delta\nu=1.5$ MHz, corresponding to $\pm 317$ km s*-1*).

We estimate the error of the HI mass measurement through jackknife resampling. From the total sample of $n$ spectra, $n/20$ randomly selected spectra are removed at a time to construct 20 jackknife samples, from which 20 mass spectra are obtained.

The jackknife estimate of the true variance of the measured value of the mass spectrum at a given frequency is then given by:

[TABLE]

where the $\langle{m}_{\rm HI}\rangle$ refers to the averaged HI mass spectrum from the original sample.

We can also measure $\langle M_{\rm HI}/L\rangle$ and its error by stacking the individual $M_{\rm HI}/L$ spectra. We do this via Equation 5 and 6, with $M_{\rm HI}$ replaced by $M_{\rm HI}/L$.

3.2 Weighting

In order to investigate the effect of different weights on the results, we explore the range $0\leq\gamma\leq 4$. As shown in Figure 5 and Table 1, highly significant values for $\langle M_{\rm HI}\rangle$ are obtained for all weighting parameters. $\langle M_{\rm HI}\rangle$ monotonically decreases as $\gamma$ increases, reflecting the lower HI mass of nearby galaxies. Similarly, $\langle M_{\rm HI}/L_{r}\rangle$ increases with $\gamma$, although the variation is somewhat less significant. The highest overall S/N occurs when $\gamma=1$. As shown in Table 1, the weighted mean redshift decreases with increasing $\gamma$. The measurements at $\gamma=1$ are more representative of the entire sample: larger $\gamma$ gives more weight to nearby galaxies; smaller $\gamma$ gives too much weight to low S/N ratio measurements of distant galaxies.

3.3 Aperture Size

With our relatively small synthesized beam area, many SDSS galaxies will be resolved or partially resolved in HI. The extraction radius therefore needs to be carefully chosen. Too small a radius may miss HI flux, while too large radius will unnecessarily introduce extra noise, and increase confusion from nearby galaxies. Based on determining the maximum radius prior to confusion becoming a problem (see Figure 6), we have chosen an aperture radius of 35 kpc, similar to the 30 kpc box size used by Geréb et al. (2015), whose observations had a somewhat smaller $\sim 10\arcsec$ synthesized beam.

Figure 6 shows that, for radii $<35$ kpc, the number of confused galaxies within (a) the aperture or within the synthesized beam, and (b) within 3 MHz (630 km s*-1*) remains in the range 120 – 130. However, at larger apertures, confusion increases rapidly, approximately doubling by 80 kpc. The luminosity distribution of the confused galaxies is shown in Figure 6.

The corresponding stacked values for $\langle M_{\rm HI}\rangle$ and $\langle M_{\rm HI}/L\rangle$ are shown in Figure 7. $\langle M_{\rm HI}\rangle$ increases monotonically, reflecting the finite size of the galaxy HI disks at small apertures, and the effect of confusion at large apertures. Between 35 and 80 kpc, $\langle M_{\rm HI}\rangle$ increases by 40 per cent. $\langle M_{\rm HI}/L\rangle$ is less sensitive to aperture. Values for both are given in Table 2, and show that S/N ratio for $\langle M_{\rm HI}/L\rangle$ is maximized when the aperture radius is 35 kpc.

3.4 Confusion correction

As shown above, $\sim 7$ per cent of our sample is potentially confused with neighbouring galaxies, both catalogued and uncatalogued. Although the WSRT synthesized beam is almost an order of magnitude smaller than the Arecibo beam and two orders of magnitude smaller than the Parkes beam, we can nevertheless estimate the corresponding correction factors for $\langle M_{\rm HI}\rangle$ and $\langle M_{\rm HI}/L\rangle$.

We have therefore carried out a mock stacking experiment on the TAIPAN+WALLABY simulations of da Cunha et al. (2017), who employ the state-of-the-art theoretical galaxy formation model GALFORM (Cole et al., 2000), in the version presented by Lagos et al. (2012). The latter follows the cosmic evolution of galaxies using a self-consistent two-phase interstellar medium model, in which stars form from the molecular gas content of galaxies. This model provides a physical distinction between atomic and molecular hydrogen in galaxies, and thus it is capable of predicting the evolution of these two components separately. The specific lightcones used here were produced using the $N$-body cold dark matter cosmological Millennium I (Springel et al., 2005) and II (Boylan-Kolchin et al., 2009) simulations, which in combination allow us to have a complete census of the HI masses of galaxies from the most HI-massive galaxies, down to an HI mass of $\approx 10^{6}$ M*⊙*. Two sets of lightcones were created and presented in da Cunha et al. (2017), one mimicking the selection function of TAIPAN and another one mimicking the selection function of WALLABY, with the primary aim of assessing the overlap population between the two surveys. Here, we use only the WALLABY111This is the Extragalactic All Sky HI Survey being carried out with the Australian Square Kilometer Array Pathfinder (Johnston et al., 2008). lightcones.

We extract 100 strips each of $2\aas@@fstack{\circ}5\times 40^{\circ}$, and in each strip we produce 35 pointings of radii 0$\aas@@fstack{\circ}$5. We select the galaxies located in these 35 pointings from $z=0.01\sim 0.11$ and $r\leq 17.7$. We also produce a volume-limited sub-sample as previously described. We use the same method as above to measure the $\langle M_{\rm HI}\rangle$ and $\langle M_{\rm HI}/L\rangle$, after locating the confused galaxies. We carry out the stacking using three methods:

Assuming that there is no confusion (i.e. stack the HI in the selected galaxies only); 2. 2.

Combine the HI in the sample galaxies with that of any companions with $r<17.77$; 3. 3.

Combine the HI in the sample galaxies with that of all companions.

We follow the method in Fabello et al. (2012) to model the confusion, estimating the total signal $S_{i}$ as the sum of the sample galaxy $S_{s}$ and the companions ($S_{c}$) weighted with two factors:

[TABLE]

where the $f_{1}$ and $f_{2}$ model the overlap between the sample galaxy and its companion in angular and redshift space.

The results are shown in Table 3. For the magnitude-limited sample, the value of $\langle M_{\rm HI}\rangle$ derived from stacking confused sample galaxies with $r\leq 17.77$ and stacking with all confused galaxies, are 1.3 $\pm$ 0.6 and 2.1 $\pm$ 0.7 per cent larger than the ‘correct’ result, respectively. For $\langle M_{\rm HI}/L\rangle$, the increase is 1.4 $\pm$ 0.6 and 1.7 $\pm$ 0.6 per cent, respectively. The increments for the volume-limited sub-samples are slightly more. For the real data, we will later utilise the confusion-included sample and use correction factors based on the ratios of method (i) and (iii) above, with 35kpc resolution.

3.5 PSF effects

In interferometric observations, the original HI sky is convolved with the point spread function (PSF) of the telescope. The PSF is then normally removed using a deconvolution algorithm. However, such a procedure is not possible when individual galaxy signals are below the noise level. Our stacks are therefore stacks of ‘dirty’ maps. To explore the effect of this, we again employ a simulation.

We convolve a simulated HI sky with the average PSF of WSRT. The simulated sky is based on the mock catalogue of Duffy et al. (2012), in which the Theoretical Astrophysical Observatory (TAO) was used to generate a light-cone catalogue from the semi-analytic models of Croton et al. (2006). Cold gas masses in this simulation were scaled by Duffy et al. to match the local HI mass function measured by ALFALFA (Martin et al., 2010) to ensure a realistic modelling of the local HI Universe. Galaxies with HI masses $M_{\rm HI}>10^{8.5}$ or apparent magnitudes $m_{r}<19.8$ are populated into the synthetic sky using the GALMOD routine from GIPSY.

In Figure 8, we illustrate the convolution process. The left panel is the PSF of WSRT, the central panel shows a 3 MHz slice of the simulated sky at $z=0.055$, and the right panel is the same slice after the convolution with the PSF. We can see clearly see the effect of sidelobes on the surrounding sky. To quantify this effect, we apply the same stacking method to the simulated sky and the convolved sky. We stack the spectra from 2727 galaxies located in the range $0.01<z<0.11$ with apparent $r$-band magnitudes brighter than 17.7.

In Table 4, we show the results of stacking with the original catalogue, the simulated sky and the convolved sky. For the latter two, we use an aperture with a radius of 35 kpc to extract the spectra. Directly stacking the HI mass given by the catalogue results in an averaged HI mass of 3.013 $\times 10^{9}h_{70}^{-2}$ M*⊙. Stacking the spectra of the selected galaxies in the simulated sky gives 3.021 $\times 10^{9}h_{70}^{-2}$ M⊙, higher due to confusion. Stacking the spectra obtained from the convolved sky gives 2.962 $\times 10^{9}h_{70}^{-2}$ M⊙*, lower due to the inclusion of negative sidelobes. Convolution makes the averaged integrated flux smaller by $1.7\%$, meaning that sidelobes only result in a small underestimate of the true signal.

3.6 Cosmic Variance

The Universe is only homogeneous on scales $>>100$ Mpc (Scrimgeour et al., 2012). Therefore observations in smaller regions can be affected by small-scale inhomogeneity, or cosmic variance. To assess the effect on our results, we assume the WSRT pointings are conical and we assign the ‘beam edge’ as the radius at which the normalized primary beam response equals to 0.1. At the median redshift of 0.066, the radius of this beam $r_{z=0.066,f=0.1}$ = 0.5195 deg, corresponding to 2368 kpc. This corresponds to a comoving volume of 6642 Mpc3 per pointing with the small volume at $z<0.01$ removed. The weighted noise-equivalent volume (square primary beam weighting) for each beam is 1545 Mpc3. The number of SDSS galaxies with spectroscopic redshifts in each pointing varies between 18 and 146 (see Table 5). Combining the 35 pointings together, the weighted sampled volume is $5.4\times 10^{4}$ Mpc3, which can be compared with the sampled volumes of HIPASS ($2.37\times 10^{5}$ Mpc3, Zwaan et al. (2005)) and the 100 per cent ALFALFA source catalog ($10.15\times 10^{5}$ Mpc3, Jones et al. (2018)).

A simple quantifiable measure of the cosmic variance can be obtained by examining the variance of galaxy counts in the TAIPAN+WALLABY simulation. We define $\xi[\%]=100\times\sigma_{\rm cv}/\langle N\rangle$, where the variance $\sigma^{2}_{\rm cv}=\Sigma_{i}(\langle N\rangle-N_{i})^{2}/n$, $\langle N\rangle$ is the mean galaxy count in the selected volumes, $N_{i}$ the number of galaxies in the volume $i$ and $n$ the total number of selected volumes. We randomly select 1000 strips of the same size as the WSRT strip and with the same redshift region from the simulation. In each strip we produce 35 pointings whose distributions are same as the WSRT observations. For galaxies within 0$\aas@@fstack{\circ}$5 of one of the pointing centres, we find $\xi=9.1\pm 0.3$ per cent.

For SDSS in the main region, the mean weighted number of galaxies at Declinations near $14^{\circ}$ across a similar 35 simulated pointings is $465$ (reduced from $1485$ by primary beam weighting), with a similar cosmic variance of 12%. However, the weighted number of galaxies in our sample is substantially higher at 519 (reduced from 1895 by weighting). This implies that the region observed is overdense by more than the variation expected from cosmic variance. Nevertheless, the cosmic variance across a wide field of view is much lower compared with a deep single pointing. Furthermore, normalization using the SDSS luminosity function removes first-order changes to the HI density associated with optical overdensities. However, second-order environmental effects may influence the final result.

4 Results

4.1 Individual Pointings

The magnitude-limited sample has a mean redshift of $\langle z\rangle=0.066$. The stacking results for each individual pointing are given in Table 5. Because of fewer galaxies and a smaller effective volume, the errors (estimated with jackknife method) are larger compared with the results from stacking the total sample. For the stacked mass spectra, only one stack (pointing 17) does not show a detection, three (pointings 12, 29 and 35) have unclear detections, while the remaining 30 pointings all result in clear detections. We show the stacked mass spectra in Appendix A.

4.2 All Galaxies

Stacking all mass spectra from our magnitude-limited sample results in a strong 67$\sigma$ detection, where the noise level is estimated from the jackknife sampling. We measured the HI mass of the stack in the manner described in Section 3. Integrating the spectral line over the rest frequency range of $\nu=1420.406\pm 1.5$ MHz and applying the confusion correction results in a mean HI mass $\langle M_{\rm HI}\rangle=(2.29\pm 0.13)\times 10^{9}h_{70}^{-2}$ M*⊙. The mean stacked value for the ratio $\langle M_{\rm HI}/L\rangle$ ratio results in a 56$\sigma$ detection with $\langle M_{\rm HI}/L\rangle=(0.306\pm 0.020)$ M${}_{\odot}/$L⊙. The stacked spectra are shown in Figure 9. For the volume-limited sub-sample, we obtain $\langle M_{\rm HI}\rangle=(0.844\pm 0.129)\times 10^{9}h_{70}^{-2}$ M⊙* and $\langle M_{\rm HI}/L\rangle=(0.369\pm 0.095)$ M${}_{\odot}/$L*⊙*.

4.3 Redshift Bins

The large redshift region and selection effects results in the sample properties changing with redshift. We split the sample into five redshift bins. The mean redshift of each bin is $\langle z\rangle=$ 0.024, 0.041, 0.062, 0.080 and 0.097. The sub-samples contain 155, 439, 453, 448 and 400 galaxies, respectively. All stacks result in significant detections. The derived average HI masses $\langle M_{\rm HI}\rangle$ and HI mass-to-light ratios $\langle M_{\rm HI}/L_{\rm r}\rangle$ are shown in Figure 10 and Table 6. The HI mass increases with redshift, and $M_{\rm HI}/L_{\rm r}$ decreases. Both results are explained by the fact that the samples are biased towards more luminous galaxies at higher redshift (see Figure 11).

5 Cosmic HI Density $\Omega_{\rm HI}$

5.1 Luminosity Bias

SDSS is a magnitude-limited sample and therefore many optically faint, but HI-rich galaxies are missed at higher redshift (Figure 4). This has an influence on our results for $M_{\rm HI}$ and $M_{\rm HI}/L$ and means that we sample different populations of galaxies at different redshifts. To account for the missed faint, but high $M_{\rm HI}/L$ ratio galaxies, we assume a power-law relation between $M_{\rm HI}/L$ and luminosity given by $M_{\rm HI}/L\sim L^{\beta}$. $\beta$ is obtained from stacking galaxies binned by their $r$-band luminosity. We show the results in Figure 11. There is a significant decrease of $M_{\rm HI}/L$ with increasing $L_{r}$. We find $\log(M_{\rm HI}/L)=(-0.587\pm 0.046)\log L+(5.246\pm 0.517)$. Since the sample is not complete in $r$-band luminosity at all redshifts, there is a selection effect in favour of low values of $M_{\rm HI}/L$ and high values of $L$ in this plot. However, only the slope of this line is relevant for the current purposes and the result appears to be similar to that derived from our the volume-limited sub-sample ($\beta=-0.662\pm 0.120$ – also shown in Figure 11). With this relation, a suitable correction for the $\langle M_{\rm HI}/L\rangle$ ratio is then given by Delhaize et al. (2013):

[TABLE]

where $\phi_{L}(L)$ is the luminosity function, $w_{i}$ is the weight of i-th galaxy and $N$ = 1895. We use $L_{\ast}$ and $\phi_{L}(L)$ given by Blanton et al. (2003), where $\phi_{L}(L)$ is a Schechter function of the form:

[TABLE]

with the following parameters: $\phi^{\ast}=5.11\times 10^{-3}h_{70}^{3}$ Mpc*-3*, $\log(L_{\ast}/L_{\odot})=10.36+\log h_{70}$ and $\alpha=-1.05$. Figure 12 shows the original and weight-corrected distribution of SDSS galaxies in $r$-band luminosity bins. The weight shifts the original distribution to lower-luminosity bins because nearby galaxies are given more weight than distant galaxies (most of which are bright). We find a correction factor of C1 = 1.38.

5.2 Stacked measurement of $\Omega_{\rm HI}$

We calculate the cosmic HI density $\rho_{\rm HI}$ from the $\langle M_{\rm HI}/L\rangle$ ratio of the stack and the luminosity density derived for SDSS galaxies. The luminosity density for $z=0.1$ in the $r$-band is given by $\rho_{L}=1.29\times 10^{8}h_{70}$ L*⊙* Mpc*-3* (Blanton et al., 2003) using 147,986 galaxies. Together with the correction factor C1, the HI density can be calculated according to:

[TABLE]

We then correct confusion according to the method described in Section 3.4. The correction for $\Omega_{\rm HI}$ is 1.7 $\pm$ 0.6 percent. Binning the galaxies into three redshift bins gives similar factors: 1.013 $\pm$ 0.006, 1.013 $\pm$ 0.006 and 1.038 $\pm$ 0.010 at mean redshifts of $\langle z\rangle$ = 0.038, 0.067 and 0.093, respectively.

After applying the above corrections for luminosity bias and confusion, we calculate a local density of $\rho_{\rm HI}=(5.46\pm 0.36)\times 10^{7}h_{70}$ M*⊙* Mpc*-3*. The error results from propagating errors in both the scaling factor and in $M_{\rm HI}/L$. To convert the local density to a cosmic HI density we divide by the $z=0$ critical density $\rho_{\rm c,0}=3H_{0}^{2}/8\pi G$ and find:

[TABLE]

Table 7 summarizes our measurement of $\Omega_{\rm HI}$ with the two correction factors consecutively applied. For the smaller volume-limited sub-sample, we find $\rm C1^{v}=1.15$, and

[TABLE]

with $\rho_{L}(z=0.024)=1.12\times 10^{8}h_{70}$ L*⊙* Mpc*-3* (given by Equation 14). The result is consistent with the magnitude-limited sample, but with larger measurement error due to the smaller sample.

We also compute $\Omega_{\rm HI}$ in different redshift bins, with the evolved $r$-band luminosity function. Using the Galaxy and Mass Assembly (GAMA) II survey, Loveday et al. (2015) found the sample to be well-fit with luminosity ($Q$) and density ($P$) evolution parameters introduced by Lin et al. (1999). The luminosity density $\rho_{L}$ can be parametrized as:

[TABLE]

with the Schechter luminosity function parameters in terms of magnitudes evolving as:

[TABLE]

[TABLE]

[TABLE]

where $P=1.0$ and $Q=1.03$ in the $r$-band. We use the results from Blanton et al. (2003) as the initial value for the Schechter parameters at $z_{0}$ = 0.1. The results in Table 8 show no measurable evolution in $\Omega_{\rm HI}$ from $z=0.038$ to $z=0.093$.

5.3 $\Omega_{\rm HI}$ in Luminosity Bins

We also measure $\Omega_{\rm HI}$ more directly in $r$-band luminosity bins using the relation:

[TABLE]

where $i$ refers to the $i$-th luminosity bin and $\phi(L)$ is the luminosity function. The $\langle M_{HI}/L\rangle_{i}$ can be obtained from Figure 11. The resultant HI density in $r$-band luminosity bins is shown in Figure 13. Using the fits to the data and summing the density in $r$-band luminosity bins from zero to infinity, we find:

[TABLE]

This is very close to the $\Omega_{\rm HI}$ derived from the stacking using the previous $\langle M/L\rangle$ bias correction. Integrating the fit in Figure 11 only in the region which has data, we have $\Omega_{\rm HI}=(2.67\pm 0.21)\times 10^{-4}h_{70}^{-1}$. If we directly sum up the data points from the stacked luminosity bins, rather than the fits, we find a value of $\Omega_{\rm HI}=(2.50\pm 0.76)\times 10^{-4}h_{70}^{-1}$. This is lower due to the HI associated with lower and higher luminosity bins than those observed.

5.4 Comparison with previous work

We show our results for $\Omega_{\rm HI}$ compared with other measurements at various redshifts in Figures 14 and 15. Each has been converted to a flat cosmology with H${}_{\circ}=70$ km s*-1* Mpc*-1* and $\rm\Omega_{m,0}$ = 0.3. $\Omega_{\rm gas}$ measurements using DLAs sometimes taken into account neutral Helium and contributions from Lyman-$\alpha$ absorbers with column densities $\log N(\rm HI)<20.3$. We convert $\Omega_{\rm gas}$ from DLAs to $\Omega_{\rm HI}$ using $\Omega_{\rm HI}=\delta_{\rm HI}\Omega_{\rm gas}/\mu$, where $\mu=1.3$ accounts for the mass of Helium and $\rm\delta_{\rm HI}=1.2$ estimates the contribution from systems below the DLA column density threshold.

As seen in Figures 14 and 15, all measurements at lower redshift ($z<0.5$) are in good agreement. But at the intermediate redshifts, measurements have large uncertainty. Our $\Omega_{\rm HI}$ measurement, marked as a red star, agrees with the measurements made at zero redshift but has a small error bar, large signal-to-noise ratio, and low systematics. It shows the usefulness of the stacking technique applied to interferometers to bridge the redshift gap between measurements using Damped Ly-$\alpha$ systems and estimates using direct 21-cm detections.

The value we measure for $\Omega_{\rm HI}$ in sub-samples at different redshifts shows no evolution, within the errors of the measurements. In combination with other results, it again suggests almost no HI gas evolution from $z\approx 0.4$ to the present, a time span of over 4 Gyr. However, combining all measurements, there remains a clear increase of $\Omega_{\rm HI}$ at higher redshift. We should note that the ‘blind’ HI 21cm surveys are measuring the ‘true’ $\Omega_{\rm HI}$ with the only assumption being that the HI 21cm emission is optically thin. On the other hand, HI stacking studies require galaxy redshifts, and are hence measuring $\Omega_{\rm HI}$ associated with galaxies detected in optical spectroscopic surveys. So high sample completeness is also required. SDSS appears to satisfy this criterion, but the under-representation of low-surface-brightness galaxies (0.1%) and close pairs $<55\arcsec$ (6%) may slightly skew the results, but this is not expected to be significant. $\Omega_{\rm HI}$ values from DLAs are similar to those from blind surveys, in that association of the gas with a galaxy is not a pre-requisite. However, there are a number of other biases such as dust obscuration, covering factor and lensing which may contribute uncertainty (Ellison et al., 2001; Jorgenson et al., 2006).

Many simulations have trouble reproducing the observed trend with redshift due the difficulty of resolving the various relevant gas phases (i.e. ionised, atomic and molecular gas, inside and outside galaxies). Recently, Davé et al. (2017), using a mid-size cosmological hydrodynamical simulation, MUFASA, found $\Omega_{\rm HI}=10^{-3.45}(1+z)^{0.74}$, which is close to the best-fit we find for the observations (Figure 14). Interestingly, previous hydrodynamical simulations have suggested that most of the HI in the Universe at $z\gtrsim 1.5-2$ is in the circumgalactic medium rather than the interstellar medium of galaxies (van de Voort et al., 2012). Using the Shark cosmological semi-analytic model of galaxy formation, Lagos et al. (2018) were able to predict the amount of atomic hydrogen contributed by the interstellar medium of galaxies to $\Omega_{\rm HI}$, across time (see Figure 14). The contribution from the interstellar medium of galaxies decreases with increasing redshift, in a trend that is the opposite to the overall increase deduced from observations.

The large impact parameters (42 kpc for ALMA J081740.86+135138.2, 18 kpc for ALMA J120110.26+211756.2 and 30 kpc for ALMAJ123055.50-113906.4) measured for the host galaxies of high-$z$ damped Lyman-alpha systems provides some support for this scenario (Neeleman et al., 2017; Neeleman et al., 2018).

It also suggests that spectral HI stacking of galaxies at redshifts beyond $z\approx 0.8$ can reveal differences between the HI content of the Universe that is accounted for in galaxies and that measured through absorption lines. Future stacking experiments at higher redshifts will therefore provide unique and stringent constraints for models of galaxy formation.

We also fit the relationship between $\Omega_{\rm HI}$ and redshift, assuming a power law relation, and find $\Omega_{\rm HI}=10^{-3.42}(1+z)^{0.68}$. A simpler linear fit to all $\Omega_{\rm HI}$ measurements, weighting all measurements according to their error, gives $\Omega_{\rm HI}(z)=0.000384+0.0002z$. The fit is shown in Figures 14. Most of the measurements are reasonably consistent with the fit, although the HI 21cm stacking result of Kanekar et al. (2016) and the HST archival study of Neeleman et al. (2016) lie below the trend.

6 Summary

In this paper we use an interferometric stacking technique to study the HI content of galaxies and confirm that there is little evolution in $\rm\Omega_{HI}$ at low redshift. Compared to previous studies, we are able to provide stronger constraints.

The data set is a 351-hr WSRT HI survey covering $\sim 7$ deg2 of the SDSS sky containing 1895 galaxies with SDSS redshifts in the range $0.01<z<0.11$. Using measurements of the mean HI mass-to-light ratio, we were able to bootstrap from the SDSS luminosity function to provide an accurate measurement of the cosmic HI gas content.

We have shown that interferometers such as WSRT offer significant advantages over single dish stacking measurements in terms of sensitivity, field-of-view and resolution which together maximize S/N ratio and minimize cosmic variance and confusion.

Over all galaxies in the sample, we find an average HI mass of $\langle M_{HI}\rangle=(2.29\pm 0.13)\times 10^{9}h_{70}^{-2}$ M*⊙* and HI mass-to-light ratio $\langle M_{HI}/L\rangle=(0.31\pm 0.02)$ M${}_{\odot}/$L*⊙. For a volume-limited sub-sample, we find $\langle M_{HI}\rangle=(0.84\pm 0.13)\times 10^{9}h_{70}^{-2}$ M⊙* and $\langle M_{HI}/L\rangle=(0.37\pm 0.09)$ M${}_{\odot}/$L*⊙*.

We derived the cosmic HI density $\Omega_{\rm HI}$ by stacking mass-to-light ratio for all galaxies. As SDSS is magnitude-limited, many optically faint but HI-rich galaxies are missing. To correct for this selection bias, we derive a weight factor which accounts for the different mass-to-light ratios of the sample compared with an unbiased selection of galaxies. We find $\rho_{\rm HI}=(5.46\pm 0.36)\times 10^{7}h_{70}$ M*⊙* Mpc*-3* and $\Omega_{\rm HI}=(4.02\pm 0.26)\times 10^{-4}h_{70}^{-1}$ at the mean redshift of $\langle z\rangle=0.066$. For a volume-limited sub-sample, we find $\Omega_{\rm HI}=(3.50\pm 0.90)\times 10^{-4}h_{70}^{-1}$ at the mean redshift of $\langle z\rangle=0.024$. We also derive the HI density from luminosity stacking and the SDSS luminosity function, finding $\Omega_{\rm HI}=4.01\times 10^{-4}h_{70}^{-1}$.

Rather than attempting to identify, then remove potentially confused targets, which has the effect of removing massive centrals and gas-rich satellites, we corrected for residual confusion using a simulation. We also explore the robustness of the result to the effect of WSRT sidelobes. For both effects, the corrections were found to be small.

Finally, we split our sample in three sub-samples with $\langle z\rangle$ = 0.038, 0.067 and 0.093 and find similar results. Our results agree well with previous $\Omega_{\rm HI}$ measurements from HI emission surveys, HI stacking and DLA surveys. Taken together, the results confirm that there seems to be little evolution in $\Omega_{\rm HI}$ at low redshift.

7 Acknowledgements

The WSRT is operated by ASTRON (Netherlands Foundation for Research in Astronomy) with support from the Netherlands Foundation for Scientific Research (NWO). This research made use of the ‘K-corrections calculator’ service available at http://kcor.sai.msu.ru/. We acknowledge the use of Miriad software in our data analysis (http://www.atnf.csiro.au/computing/software/miriad/). This research made use of the Sloan Digital Sky Survey archive. The full acknowledgment can be found at http://www.sdss.org. Parts of this research were supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013.

Appendix A Stacked Spectra

We show the stacked mass spectra for pointings 1 - 35 in Figure 16, Figure 17 and Figure 18. The red-dashed lines show the region over which we do the integration to compute the average HI mass. For the stacked mass spectra, only one stack (pointing 17) does not show a detection, three (pointings 12, 29 and 35) show unclear detections, while the remaining 30 pointings all result in clear detections.