The neutral hydrogen cosmological mass density at z = 5

Abstract

1 INTRODUCTION

The neutral hydrogen mass density of the Universe, |$\Omega_{\rm H\,{\small i}}$|⁠, is an important cosmological observable. It determines the precision with which cosmological parameters can be constrained by observations of the H i intensity power spectrum (e.g. Barkana & Loeb 2007; Chang et al. 2008; Wyithe & Loeb 2008; Padmanabhan, Choudhury & Refregier 2015), and we expect its evolution to be linked to the cosmic star formation history. The main contributor to |$\Omega_{\rm H\,{\small i}}$| is high column density, predominantly neutral gas clouds (e.g. O'Meara et al. 2007; Zafar et al. 2013), self-shielded from ionizing radiation and therefore likely fuel for future star formation (e.g. Wolfe, Gawiser & Prochaska 2005). Thus, tracing the evolution of |$\Omega_{\rm H\,{\small i}}$| from the end of reionization, through the epoch of the cosmic star formation peak at z ∼ 2 to the present day, is of central importance to our understanding of galaxy formation. It also provides an excellent integral constraint against which theoretical models of galaxy formation can be tested.

At redshift <0.3, H i 21 cm emission can be used to measure |$\Omega_{\rm H\,{\small i}}$| either directly or by stacking analyses (e.g. Zwaan et al. 2005; Martin et al. 2010). At higher redshifts, where emission is too weak to be detected with current facilities, |$\Omega_{\rm H\,{\small i}}$| can instead be inferred from the incidence rate of damped Lyα systems (DLAs, defined as absorption systems with |${\it N}_{\rm H\,{\small i}}\ge 20.3$| cm⁻²), which trace the bulk of neutral gas in the Universe (Prochaska, Herbert-Fort & Wolfe 2005). These systems are detected in absorption in the spectra of background QSOs, and their characteristic damping wings allow column densities to be measured even at low spectral resolution.

Early DLA surveys at 2 < z < 4, which were typically comprised of a few hundred QSOs and assumed a cosmological deceleration parameter q₀ = 0.5 or 0, suggested that the gas mass density in DLAs may have been sufficient to produce most of the stars seen in the local Universe (Lanzetta et al. 1991; Wolfe et al. 1995; Storrie-Lombardi, McMahon & Irwin 1996). However, a change to a modern concordance cosmology revealed that DLAs at z ∼ 3 contain <50 per cent of the present-day mass density in stars (e.g. Storrie-Lombardi & Wolfe 2000; Péroux et al. 2005, see also Section 5.2). In addition, recent DLA surveys at 2 < z < 4 using more than 10 000 QSOs assembled from the Sloan Digital Sky Survey (SDSS; Prochaska & Herbert-Fort 2004; Prochaska et al. 2005; Noterdaeme et al. 2009, 2012; Prochaska & Wolfe 2009) have shown that there is very little evolution in the H i mass density from z = 3 to the present day. This is starkly at odds with the strong evolution in the star formation rate over the same period (e.g. Madau & Dickinson 2014). One view is that H i represents a transitory phase fuelling star formation (e.g. Prochaska et al. 2005; Davé et al. 2013), which is continually replenished by more highly ionized gas from either the intergalactic medium (IGM) or recycled galactic outflows.

While it is important to constrain |$\Omega_{\rm H\,{\small i}}$| across the whole of cosmic history, it is of particular interest at the highest redshifts. Rafelski et al. (2014) report a decrease in the metal mass density in DLAs from z = 5 to 4.5, hinting at an abrupt change in the enrichment of H i gas past z = 5. This may be caused by a change in the population of objects containing neutral hydrogen, which could be accompanied by a similarly abrupt evolution in |$\Omega_{\rm H\,{\small i}}$|⁠. Moreover, since massive stars in galaxies are believed to have reionized the Universe (e.g. Bouwens et al. 2012), it is important to track the evolution of the fuel for star formation up to the epoch of reionization. However, it is a challenge to assemble the large sample of high-redshift QSO spectra necessary for a z > 4.5 DLA survey. The decline in the QSO space density at z > 3 means that relatively few redshift >4.4 QSOs were observed by the SDSS, and those that were typically have too low a signal-to-noise (S/N) to reliably identify DLAs. For example, Rafelski et al. (2012, 2014) find a misidentification rate of 26 per cent for DLA candidates from SDSS DR5 at z > 4 and of 97 per cent for candidates from DR9 at z > 4.7. For this reason, smaller DLA surveys have been performed at higher redshift, often using higher resolution spectra to make robust identifications of DLAs. Péroux et al. (2003), Guimarães et al. (2009) and Songaila & Cowie (2010) have all presented measurements of |$\Omega_{\rm H\,{\small i}}$| at z > 4.5. Songaila & Cowie (2010, hereafter S10) give a cumulative result including data from all these previous studies, and this represents the highest redshift measurement of |$\Omega_{\rm H\,{\small i}}$| to date. They use a sample of 19 QSOs with emission redshifts >4.5, and their measurement hints at a possible downturn in |$\Omega_{\rm H\,{\small i}}$| at z ≥ 4, but the uncertainties from sample variance at z > 4.3 are large.

Here we measure |$\Omega_{\rm H\,{\small i}}$| as traced by DLAs at 3.5 < z < 5.4 using a homogeneous sample of 163 QSOs with emission redshifts between 4.4 and 5.4. This represents an increase in redshift path of a factor of 8 over S10 at z > 4.5. Identifying DLAs becomes increasingly difficult at higher redshift, as H i absorption from the highly ionized IGM becomes more severe, and blending with strong systems below the DLA threshold can cause misidentification of DLAs. Therefore, we carefully check for systematic misidentifications in our sample using both mock spectra and higher resolution spectra of DLA candidates. More than 70 per cent of our DLA candidates (and >85 per cent at z > 4.5) have been observed at higher resolution (Rafelski et al. 2012, 2014), allowing us to confirm their |${\it N}_{\rm H\,{\small i}}$| despite the increased IGM blending at high redshift.

This paper is structured as follows. In Section 2, we describe the QSO spectra used for the analysis. Section 3 describes the formalism used to derive |$\Omega_{\rm H\,{\small i}}$| from our observations, and Section 4 describes our method for measuring the DLA incidence rate, accounting for systematic effects. Section 5 describes our main result, a measurement of the neutral hydrogen mass density at z = 5, and discusses its implications. Section 6 summarizes our conclusions. We assume a flat Λ cold dark matter cosmology, with H₀ = 70 km s^{− 1} Mpc⁻¹, Ω_{m, 0} = 0.3 and Ω_{Λ, 0} = 0.7. All distances are comoving unless stated otherwise. The data and code used for this paper are available at https://github.com/nhmc/GGG_DLA.

2 DATA

Figure 1.

Open in new tabDownload slide

The left-hand panel shows the emission redshift distribution for QSOs in the low-resolution GGG sample (open histogram), and for the subsample of these QSOs targeted with higher resolution spectra. The right-hand panel shows the redshift path, g(z), for detecting DLAs for the GGG sample. g(z) is defined as the number of QSOs where a DLA can be detected as a function of DLA redshift. For comparison, g(z) for previous high-redshift DLA surveys are also shown: for Péroux et al. (2003) and for z > 4.5 QSOs from S10. We do not show g(z) for the SDSS DLA surveys (e.g. Noterdaeme et al. 2012). Their g(z) formally extends to z > 4, but Prochaska et al. (2005) warn that this high-redshift sensitivity should be viewed conservatively and Noterdaeme et al. (2012) do not include DLAs with z > 3.5 in their statistical sample.

For a detailed description of the GGG spectra and the procedure used to reduce them, see Worseck et al. (2014). In brief, they were observed with the Gemini Multi-Object Spectrometers on the Gemini telescopes, yielding a typical S/N ∼ 20 per 1.85 Å pixel in the Lyα forest at a resolution of ∼5.5 Å (full width at half-maximum, FWHM). The spectral coverage was tuned to be roughly constant in the quasar rest frame (typically 850–1450 Å). The high-resolution ESI spectra we use¹ have a typical S/N of 15 per 10 km s^{− 1} pixel and a resolution FWHM of 31 km s^{− 1} (see Table 1). The single MagE spectrum has a similar S/N but a resolution of 56 km s^{− 1}.

3 FORMALISM

Our aim is to measure the cosmic H i mass density at 3.5 < z < 5.4. The bulk of the neutral gas at 2 < z < 5 is in DLAs, with an ∼15 per cent contribution from sub-DLAs (which have |$10^{19}<{\it N}_{\rm H\,{\small i}}/(\rm cm^{-2})<10^{20.3}$|⁠) and more highly ionized Lyman limit and Lyα forest absorbers with |${\it N}_{\rm H\,{\small i}}<10^{19}\,\rm cm^{-2}$| (Péroux et al. 2005; Prochaska et al. 2005; O'Meara et al. 2007; Zafar et al. 2013). There are several ways to express the comoving mass density of neutral hydrogen used in the literature. For measurements at low redshift using radio emission, authors typically quote |$\Omega _{\rm H\,{\small i}}$|⁠, which is the mass of neutral hydrogen alone, excluding any mass in molecules and helium. For DLA absorption studies, authors generally quote the gas mass in DLAs, |$\Omega ^\mathrm{DLA}_{g}$| (sometimes the g subscript is omitted) including a factor μ to account for helium. Prochaska et al. (2005) advocate using the quantity |$\Omega ^\mathrm{neut}_g$|⁠, which is the mass in predominantly neutral gas, which can be different from |$\Omega ^\mathrm{DLA}_{g}$|⁠. In this work, we quote the mass density from H i alone, |$\Omega_{\rm H\,{\small i}}$|⁠, and exclude any mass contribution from helium or molecules. Due to contamination and the low resolution of the GMOS spectra, we only measure H i in DLAs, |$\rm \Omega ^{\rm DLA}_{\rm H\,\small {I}}$|⁠. To convert to |$\Omega_{\rm H\,{\small i}}$|⁠, we apply a correction derived from measurements of lower |${\it N}_{\rm H\,{\small i}}$| systems in previous work.

We measure |$\Omega^{\rm DLA}_{\rm H\,{\small i}}$| by counting the incidence rate of DLAs in the spectra and measuring |${\it N}_{\rm H\,{\small i}}$| from their strong damping wings. Below is a summary of the formalism used to derive |$\Omega ^{\rm DLA}_{\rm H\,\small {I}}$| from the DLA incidence rate. See section 4.1 of Prochaska et al. (2005) and the review by Wolfe et al. (2005) for a more detailed description.

3.1 Other systematic effects contributing to |$k({\it N}_{\rm H\,{\small i}})$|

In measuring |$k({\it N}_{\rm H\,{\small i}})$|⁠, we explicitly take into account the rate of spurious DLAs (false positives) and missed DLAs (false negatives). There are several other systematic effects which could also contribute to |$k({\it N}_{\rm H\,{\small i}})$|⁠, which we discuss here.

The first of these is any uncertainty in the |${\it N}_{\rm H\,{\small i}}$| measurements. If there are large uncertainties in |${\it N}_{\rm H\,{\small i}}$|⁠, or systematic offsets in the |${\it N}_{\rm H\,{\small i}}$| estimated from the spectra as a function of |${\it N}_{\rm H\,{\small i}}$|⁠, this may change the inferred |$f({\it N}_{\rm H\,{\small i}})$|⁠. However, in Section 4.3 we show that the |${\it N}_{\rm H\,{\small i}}$| error from the GMOS spectra (0.2 dex) does not have a detectable systematic bias, and Section 5 shows that any errors it introduces to |$\Omega_{\rm H\,{\small i}}$| are negligible compared to other uncertainties. A related effect is for |${\it N}_{\rm H\,{\small i}}$| measurements at the DLA threshold of |${\it N}_{\rm H\,{\small i}}$|=10^20.3 cm⁻², where the more numerous lower column density systems may be counted as DLAs through |${\it N}_{\rm H\,{\small i}}$| uncertainties. This bias is a net source of false positives, and so should be taken into account by our procedure for estimating k_real.

A second possibility is the presence of dust in DLAs. If DLAs contain large amounts of dust, they are able to extinguish the light from a background QSO, removing these sightlines from our survey. In this case, we would measure a lower incidence of high-metallicity, high-|${\it N}_{\rm H\,{\small i}}$| DLAs, which presumably contain the most dust. However, several studies have shown that most DLAs are not associated with significant amounts of dust (e.g. Murphy & Liske 2004; Vladilo, Prochaska & Wolfe 2008), and DLAs towards radio-selected QSOs, which are insensitive to the presence of dust, have a similar |${\it N}_{\rm H\,{\small i}}$| distribution to those in optically selected QSOs (Ellison et al. 2001; Jorgenson et al. 2006). Pontzen & Pettini (2009) find that the cosmic H i mass density may be underestimated by 3–23 per cent at z ∼ 3 due to selection biases from dust. We do not include this relatively small effect in our analysis, but note where its inclusion would affect our conclusions.

Gravitational lensing may also introduce a bias. DLA host galaxies may lens background QSOs, making them more likely to be found in our survey. This would result in brighter QSOs being more likely to show foreground DLA absorption compared to fainter QSOs. At z ∼ 3, Murphy & Liske (2004) found evidence at the ∼2σ level that DLAs tend to be found towards brighter QSOs. Prochaska et al. (2005) found a higher incidence rate of high-|${\it N}_{\rm H\,{\small i}}$| DLAs towards brighter QSOs compared to fainter QSOs over a redshift range 2–4.5, which resulted in a significant (>95 per cent) difference in |$\Omega_{\rm H\,{\small i}}$| between the two samples. They attributed this effect to gravitational lensing. We confirm that this effect is also present in our sample (which has some overlap with the Prochaska et al. sample): there is a 25 ± 15 per cent higher incidence rate of DLAs towards QSOs with z-band magnitude ≤19.2 compared to QSOs with z > 19.2 mag. DLAs towards bright QSOs also tend to have high |${\it N}_{\rm H\,{\small i}}$|⁠, resulting in a 30 per cent increase in |$\Omega_{\rm H\,{\small i}}$| for the brighter compared to the fainter QSO sample. The significance of the excesses we measure is modest (1.7σ), and a Kolmogorov–Smirnov test between the |${\it N}_{\rm H\,{\small i}}$| distributions towards z ≤ 19.2 and z > 19.2 mag quasars yields D = 0.3 and a probability of 22 per cent that the two samples are drawn from the same underlying distribution. Therefore, while this difference hints at a selection effect related to the background QSO brightness, we cannot yet rule out a simple statistical fluctuation. We further discuss how this possible bias may affect our |$\Omega_{\rm H\,{\small i}}$| measurement in Section 5.1.1.

3.2 Conversion from |$\Omega ^{\rm DLA}_{\rm H\,\small {I}}$| to |$\Omega _{\rm H\,\small {I}}$|

Previous absorption studies have shown that the dominant contribution to |$\Omega_{\rm H\,{\small i}}$| is from DLAs. Lower column density systems also contribute an appreciable fraction of |$\Omega_{\rm H\,{\small i}}$|⁠, however. This fraction is 15–30 per cent at z = 3, depending on the assumed |${\it N}_{\rm H\,{\small i}}$| distribution (e.g. O'Meara et al. 2007; Noterdaeme et al. 2009; Prochaska, O'Meara & Worseck 2010; Zafar et al. 2013). To parametrize this uncertainty, we introduce a correction factor |$\delta _{\rm H\,\small {I}}\equiv \rm \Omega _{\rm H\,\small {I}}/\rm \Omega ^{\rm DLA}_{\rm H\,\small {I}}$| to convert between |$\rm \Omega ^{\rm DLA}_{\rm H\,\small {I}}$|⁠, which we measure, and |$\Omega_{\rm H\,{\small i}}$|⁠. We assume that the |${\it N}_{\rm H\,{\small i}}$| distribution at z > 4 is not dramatically different from that at z ∼ 3 and take δ = 1.2, which implies a 20 per cent contribution from lower column density systems. Zafar et al. find that the contribution of sub-damped systems to |$\Omega_{\rm H\,{\small i}}$| increases with redshift, possibly due to a weakening of the UV background as the number density of QSOs drops at high redshift. Therefore, a goal of future surveys should be to measure the contribution of these sub-damped systems at z > 4.

4 METHOD

4.1 Procedure for identifying DLAs

We measure the frequency of DLAs, f_meas, by identifying DLA candidates by eye in the GMOS spectra, and then correcting for any biases in identification using mock spectra. To identify candidates, we performed the following steps for each QSO spectrum.

• Estimate the continuum as a spline, placing the spline knot points by hand. We used the low-z composite QSO spectrum from Shull, Stevans & Danforth (2012) to indicate the position of likely QSO emission lines which fall inside the Lyα forest.

• Look for a possible damped Lyα line in the Lyα forest between the QSO Lyα and Lyβ emission lines. Estimate its redshift and |${\it N}_{\rm H\,{\small i}}$| by plotting a single-component Voigt profile with b = 30 km s^{− 1} over the spectrum, and varying |${\it N}_{\rm H\,{\small i}}$| and z until it matches the data by eye.² If necessary the continuum was varied at the same time |${\it N}_{\rm H\,{\small i}}$| was estimated to obtain a plausible fit. At higher redshifts, blending with IGM absorption can make estimating |${\it N}_{\rm H\,{\small i}}$| challenging, as the damping wings can be very heavily blended with IGM absorption. In this case, the best constraint on |${\it N}_{\rm H\,{\small i}}$| is not from the shape of the damping wings, but instead from the extent of the Lyα trough consistent with zero flux and from any higher order Lyman transitions.

• If a candidate DLA is found based on the Lyα profile, use its higher order Lyman series (if available in the spectrum) to refine its redshift and |${\it N}_{\rm H\,{\small i}}$|⁠.

• Repeat steps (ii) and (iii) for all DLA candidates in the Lyα forest.

DLA absorption can also be detected bluewards of the QSO Lyβ emission line. However, we chose to search only between Lyα and Lyβ emission in our sample to maximize the chance of having useful Lyman series lines in addition to Lyα, and to avoid any additional systematic effects caused by further blending with the Lyβ forest. While most DLAs also have associated metal lines detected by the GMOS spectra, we did not use any metal line information when measuring the DLA candidate redshift or |${\it N}_{\rm H\,{\small i}}$|⁠. This was done to avoid any bias against finding low-metallicity systems, which may not have detectable metals in the GMOS spectra.

Figure 2.

Open in new tabDownload slide

DLAs identified in the GMOS spectra (resolution FWHM ∼230 km s^{− 1}) which are confirmed in higher resolution ESI spectra (resolution FWHM ∼30 km s^{− 1}). In each case, the top panels show the GMOS spectrum and the bottom panels the ESI spectrum of the same QSO. The model shows the |${\it N}_{\rm H\,{\small i}}$| and redshift estimated from the ESI spectra with the redshift fixed by low-ion metal absorption. The shaded region shows an uncertainty in |$\log {\it N}_{\rm H\,{\small i}}$| of 0.2. The |${\it N}_{\rm H\,{\small i}}$| and redshift estimated from the GMOS spectra are given in Table 2.

With these DLA candidate lists, we can derive the measured incidence rate of DLAs, f_meas. However, despite our attempt to take continuum uncertainties and IGM absorption into account when measuring |${\it N}_{\rm H\,{\small i}}$| for each DLA, large systematic uncertainties may remain. The following sections describe how we quantify these uncertainties using the correction factors k_real and k_found to f_meas.

4.2 Estimation of k_real and k_found

We expect k_real to be less than unity, meaning that there are some spurious DLA candidates. The rate of these spurious candidates is estimated in two ways. First, we use the sample of higher resolution spectra to identify DLAs, and compare these with the DLA candidates found in the low-resolution sample. Secondly, we create mock low-resolution spectra which closely match the GMOS spectra and contain DLAs generated from a distribution at z = 3, and then search these spectra for DLAs in the same way as the real spectra.

k_found is also expected to be less than unity, which means some true DLAs exist which we do not flag as DLA candidates in the low-resolution spectra. Again we estimate the fraction of true DLAs recovered in two independent ways, using higher resolution spectra and mocks. In the first case, DLAs identified in the higher resolution QSO spectra were used as a reference list of true DLAs, and compared to the candidate DLAs found in the lower resolution spectra of the same QSOs. In the second case, we used mock GMOS spectra, which allow us to directly compare known DLAs in the spectra to the DLA candidates.

Our motivation for using two different ways to estimate the correction factors (mocks and high-resolution spectra) is to test different systematic effects. The main advantage of the mocks is that the true DLA properties are known precisely. However, while we attempt to reproduce the real spectra as closely as possible, including Lyα forest clustering, QSO redshift and S/N distribution, it is still possible that the mocks may differ from the real GMOS spectra. Metal absorption (not included in the mocks) or clustering of strong absorbers that is different from the mocks may cause more spurious DLAs. Alternatively, non-Gaussian noise in the real spectra at low fluxes may mean that true DLAs are more likely to be missed in the real spectra. Conversely, for the high-resolution sample, the true DLA properties are not known with complete certainty, but the correct clustering, IGM blending, noise and metal absorption are all included. Therefore, these two approaches provide complementary estimates of k_found and k_real. The following sections describe these approaches in more detail.

4.2.1 Corrections using high-resolution spectra

DLAs can be found more easily in our sample of high-resolution spectra, and their |${\it N}_{\rm H\,{\small i}}$| and redshift are more accurately measured, in comparison to the lower resolution GMOS spectra. Therefore, we independently identify DLAs in these spectra for the purpose of deriving the correction factors k_real and k_found, and to test for any systematics in estimating |${\it N}_{\rm H\,{\small i}}$| and z for each DLA. When identifying the DLAs in the 59 high-resolution spectra, we follow the same process outlined for the lower resolution spectra in Section 4.1, using the Lyman series to estimate the redshift and |${\it N}_{\rm H\,{\small i}}$|⁠. However, we also refine the redshift and |${\it N}_{\rm H\,{\small i}}$| using the position of low-ionization metal lines (O i, Si ii, C ii and Al ii) where possible. For the 20 QSOs with high-resolution spectra which are not in the GGG sample, we created low-resolution spectra by convolving the high-resolution spectra to the same FWHM resolution, and rebinning to the same pixel size as the GMOS spectra. The same noise array was used for these spectra as for the GGG QSO with a redshift closest to each QSO, normalizing such that the median S/N within rest-frame wavelengths 1260–1280 Å match. These low-resolution spectra were searched for DLAs in the same way as the GMOS spectra.

In this way, we made two lists of DLAs, one from the high-resolution spectra and another from low-resolution spectra of the same QSOs. The DLAs identified in the higher resolution sample are listed in columns 5 and 6 of Table 2. We then estimated k_real as N_{cand, true}/N_cand, where N_cand is the number of DLA candidates from the low-resolution spectra and N_{cand, true} is the number of those candidates confirmed to be DLAs by the high-resolution spectra. k_found is estimated as N_{cand, true}/N_true, where N_true is the number of DLAs found in the high-resolution spectra and N_{cand, true} is the number of those also flagged as DLA candidates in the low-resolution spectra. We calculate the binomial confidence intervals on k_real and k_found using the method described by Cameron (2011).

With this procedure, we find |$k_\mathrm{real}= 0.80_{-0.08}^{+0.07}$| and |$k_\mathrm{found}=0.84^{+0.06}_{-0.08}$| using DLAs identified by JXP (see Figs 5 and 6) with similar values found by NHMC. Both are below unity, and so there are both spurious DLA candidates and real DLAs missed. Spurious DLAs usually occur when flux spikes are smoothed away at GMOS resolution, making a lower |${\it N}_{\rm H\,{\small i}}$| system appear to have strong damping wings. An example spurious DLA is shown in Fig. 3. Real DLAs are generally missed due to flux fluctuations in the core of the Lyα line: an example is shown in Fig. 4.

Figure 3.

Open in new tabDownload slide

Example of a spurious DLA candidate. This was identified as a DLA with |${\it N}_{\rm H\,{\small i}}=10^{20.4\pm 0.2}$| cm⁻² in the GMOS spectrum shown in the top panels. However, the residual flux spikes at Lyα and Lyγ in the higher resolution (FWHM ∼30 km s^{− 1}) ESI spectrum in the bottom panels show that this system must have |${\it N}_{\rm H\,{\small i}}<10^{20.3}$| cm⁻².

Figure 4.

Open in new tabDownload slide

A DLA that was not correctly identified in the GMOS spectra. Lower panels show the DLA in the ESI spectrum, with |${\it N}_{\rm H\,{\small i}}=10^{20.55\pm 0.2}$| cm⁻². The residual flux in the core of the Lyα line in the GMOS spectrum, however (top-left panel), meant this system was missed. This residual flux around velocities Δv ∼ 0 km s^{− 1} may be caused by either statistical fluctuations or systematics associated with sky background level.

4.2.2 Corrections using mock spectra

Our method for generating mock spectra is described in Appendix A. In this case, the |${\it N}_{\rm H\,{\small i}}$| for each DLA is known, and so can be directly compared to the candidates identified in the low-resolution mocks. Again k_real is estimated as N_{cand, true}/N_cand, where N_cand is the number of DLA candidates from the low-resolution mock spectra and N_{cand, true} is the number of those candidates that are DLAs. k_found is estimated as N_{cand, true}/N_true, where N_true is the true number of DLAs in the mocks and N_{cand, true} is the number of those recovered as DLA candidates. Again we calculate the errors on k_real and k_found assuming a binomial confidence interval. For the mocks, we find k_real = 0.71 ± 0.06 and |$k_\mathrm{found}=0.92^{+0.04}_{-0.07}$| using DLAs identified by JXP (see Figs 5 and 6). Similar values are found by NHMC (see Figs A2 and A3).

Figure 5.

Open in new tabDownload slide

The fraction of true DLAs that were correctly identified by one of the authors (JXP), k_found, as a function of the true redshift and |${\it N}_{\rm H\,{\small i}}$|⁠. Top panels are for the high-resolution sample, and bottom panels are for mocks. The upper row of numbers under each histogram gives the number of DLA candidates that are correct per bin, and the lower row the total number of candidates. The total numbers for all bins are given at the top in the left-hand panels. Vertical lines show the binomial 68 per cent uncertainties.

Figure 6.

Open in new tabDownload slide

The fraction of non-spurious DLA candidates, k_real, by one of the authors (JXP), as a function of the candidate redshift and |${\it N}_{\rm H\,{\small i}}$|⁠. Top panels are for the high-resolution sample, and bottom panels are for mocks. The upper row of numbers under each histogram gives the number of true DLAs that are recovered per bin, and the lower row the total number of true DLAs. The total numbers for all bins are given at the top of the left-hand panels. Vertical lines show the binomial 68 per cent uncertainties.

4.2.3 Comparison of correction factors and their dependence on redshift and column density

We expect k_real and k_found to be a function of a DLA's |${\it N}_{\rm H\,{\small i}}$| (high-|${\it N}_{\rm H\,{\small i}}$| candidates should be more reliable), spectral S/N (low-S/N spectra will produce more spurious candidates) and redshift (more spurious DLAs will be found at high redshift where there is more IGM absorption). The most important of these for our measurement of |$\Omega_{\rm H\,{\small i}}$| is any redshift or |${\it N}_{\rm H\,{\small i}}$| dependence. Noterdaeme et al. (2009, 2012) show that at z ∼ 2.5, systems with |${\it N}_{\rm H\,{\small i}}=10^{20.6{\rm -}21.5}\,\rm cm^{-2}$| make the largest contribution to |$\Omega_{\rm H\,{\small i}}$|⁠. Thus, we expect completeness corrections in this column density range to have the largest effect on the final derived |$\Omega_{\rm H\,{\small i}}$|⁠.³

The top panels of Fig. 5 show the correction factor k_found from the high-resolution spectra binned by the true DLA redshift and |${\it N}_{\rm H\,{\small i}}$|⁠, and the bottom panels show the same correction factor estimated from the mocks. Fig. 6 shows the correction factor k_real binned by the candidate DLA redshift and |${\it N}_{\rm H\,{\small i}}$|⁠, again for the high-resolution spectra and mocks. These are derived from DLAs identified by one of the authors (JXP) who search the spectra for DLAs, but values for the other author (NHMC) are similar. There is no evidence for a strong dependence of k_real or k_found on redshift, using either the high-resolution spectra or the mocks. However, there is a weak dependence of k_real and k_found on |${\it N}_{\rm H\,{\small i}}$|⁠, with the lowest |${\it N}_{\rm H\,{\small i}}$| bin having a significantly lower k_real than for higher |${\it N}_{\rm H\,{\small i}}$| bins. This matches our expectations: weaker candidate DLAs are more likely to be spurious, and true DLAs that are weak are more likely to be missed. We take this |${\it N}_{\rm H\,{\small i}}$| dependence into account when applying the correction factors as described in Section 5. We find no strong dependence of the correction factors on S/N in either the mocks or the high-resolution sample for the range of S/N the GMOS spectra cover.

Figs 5 and 6 also show that corrections derived from the mocks and high-resolution spectra are in reasonable agreement. The main difference is in the number of spurious systems with |${\it N}_{\rm H\,{\small i}}\sim 10^{20.3{\rm -}20.6}$| cm⁻². The right-hand panels of Fig. 6 show that there are more weak, spurious DLAs found in the mocks compared to the real GMOS spectra. However, we show in the following section that the correction factor in this |${\it N}_{\rm H\,{\small i}}$| range is not important for estimating |$\Omega_{\rm H\,{\small i}}$|⁠, and for the remaining bins the mocks and high-resolution corrections match to within 20 per cent. As we discussed earlier, the high-resolution sample and mocks test different systematic uncertainties which may affect |$\Omega_{\rm H\,{\small i}}$|⁠. Therefore, the consistency of the correction factors between these two methods suggests that the mocks reproduce the true GMOS spectra well, and that DLAs have been identified correctly in the higher resolution spectra.

4.3 Uncertainties in |${\it N}_{\rm H\,{\small i}}$| and redshift

If DLA column densities estimated from the GMOS spectra are systematically in error, our measurement of |$\Omega_{\rm H\,{\small i}}$| may be biased. Such a systematic could occur because of incorrect placement of the continuum or blending of damping wings with the Lyα forest. This is an additional effect not accounted for by the correction factor, k, to f_meas. Therefore, we search for any systematic offset in |${\it N}_{\rm H\,{\small i}}$| by matching DLA candidates from the low-resolution spectra to known DLAs in the high-resolution sample and mocks.

The results of this test are shown in Fig. 7. The log |${\it N}_{\rm H\,{\small i}}$| difference is plotted as a function of redshift, the true |${\it N}_{\rm H\,{\small i}}$| and S/N for the high-resolution sample (top panels) and mocks (bottom panels). For both the mocks and high-resolution samples, both |$\Delta \log {\it N}_{\rm H\,{\small i}}$| and Δv are centred on 0. The standard deviation of the velocity and |$\log {\it N}_{\rm H\,{\small i}}$| offsets are 184/216 km s^{− 1} and 0.165/0.196 for the high-resolution sample and mocks, respectively. We therefore adopt 0.2 dex as our uncertainty in |${\it N}_{\rm H\,{\small i}}$|⁠. There is no trend seen with redshift or S/N. There may be a trend with |${\it N}_{\rm H\,{\small i}}$|⁠, but above N = 10^20.3 cm^{− 2}, it is too weak to significantly affect |$\Omega_{\rm H\,{\small i}}$|⁠. We conclude that there is no systematic bias in |${\it N}_{\rm H\,{\small i}}$| which might adversely affect the |$\Omega_{\rm H\,{\small i}}$| measurement.

Figure 7.

Open in new tabDownload slide

The difference between |${\it N}_{\rm H\,{\small i}}$| estimated for DLA candidates in low-resolution spectra (⁠|$N_\mathrm{H\,\small {I}, cand}$|⁠) and |$N_\mathrm{H\,\small {I}, true}$| measured from high-resolution spectra (top) or known from mock line lists (bottom), and similarly for the velocity offset from the true DLA redshift. This is for one of the authors (JXP), but the results for NHMC are similar. These show that there is no strong systematic offset in the estimated |${\it N}_{\rm H\,{\small i}}$| as a function of redshift, S/N or |${\it N}_{\rm H\,{\small i}}$| which might systematically bias |$\Omega_{\rm H\,{\small i}}$| significantly. Grey shading shows the regions that cannot be populated due to the requirement that both |$N_\mathrm{H\,\small {I}, cand}$| and |$N_\mathrm{H\,\small {I}, true}$| be >10^20.3 cm⁻².

Fig. 7 also shows the redshift difference between matched DLAs expressed as a velocity difference. DLAs identified in the higher resolution spectra use low-ionization metal lines to set a precise DLA redshift with an error of a few km s^{− 1}. Both the mocks and high-resolution sample show that an uncertainty of ∼200 km s^{− 1} results from estimating redshifts using Lyman series absorption alone (without reference to metal absorption) in the low-resolution spectra.

4.4 DLA incidence rate and differential |${\it N}_{\rm H\,{\small i}}$| distribution

Fig. 8 shows the differential |${\it N}_{\rm H\,{\small i}}$| distribution from the GGG sample compared to that from the SDSS sample from Prochaska & Wolfe (2009), which is consistent with the more recent estimate from Noterdaeme et al. (2012). We have four different measurements of the correction factor |$k({\it N}_{\rm H\,{\small i}})$|⁠, from two different authors using the mocks and high-resolution spectra, so there are four different estimates of |$f({\it N}_{\rm H\,{\small i}},X)$|⁠. We find the final |$f({\it N}_{\rm H\,{\small i}},X)$| by averaging these four estimates. The uncertainties on this value include a statistical and systematic component. The statistical uncertainty is found by bootstrap resampling, using 1000 samples from the observed DLA distribution, and averaging these uncertainties for the four different estimates. The systematic uncertainty is then assumed to be the standard deviation in the four estimates. These systematic and statistical components are added in quadrature to give the errors shown in Fig. 8. The two distributions are similar overall, although there is a clear discrepancy between the GGG and z = 3 |$f({\it N}_{\rm H\,{\small i}},X)$| for the bin at |$\log {\it N}_{\rm H\,{\small i}}\sim 21.2$|⁠, which hints at evolution in the shape of |$f({\it N}_{\rm H\,{\small i}},X)$| at high redshift. However, a simple change in the normalization is also consistent with the data.

Figure 8.

Open in new tabDownload slide

The column density distribution, |$f({\it N}_{\rm H\,{\small i}})$|⁠, for the GGG sample of DLAs. Black points show the measurements with correction |$k({\it N}_{\rm H\,{\small i}})$| applied. Both are slightly offset in |${\it N}_{\rm H\,{\small i}}$| for clarity. Red points show the SDSS DR5 measurements from Prochaska & Wolfe (2009) after applying the correction to the redshift search path recommended by Noterdaeme et al. (2009). The errors are 1σ, and include statistical and systematic errors (see Section 4.4 for more details).

The DLA incidence rate, ℓ(X), is shown in Fig. 9. This observable is more sensitive to the lowest |${\it N}_{\rm H\,{\small i}}$| DLAs than |$\Omega_{\rm H\,{\small i}}$|⁠. Since the correction factors we derive are strongest for low-|${\it N}_{\rm H\,{\small i}}$| DLAs and these DLAs have a strong effect on ℓ(X), we expect ℓ(X) to be sensitive to the particular choices of correction factors. This is indeed the case – there are systematic differences at least as large as the statistical errors, and they depend on whether the mocks or the high-resolution spectra are used to estimate the correction factor. Similarly, large differences are found between ℓ(X) by each of the two authors who searched for DLAs. The ℓ(X) values we measure are consistent with a smooth increase from z = 2 to 5. However, since we do not know which |$k({\it N}_{\rm H\,{\small i}})$| correction factors are best, we do not attempt to present a definitive ℓ(X) measurement here. A large sample of higher resolution spectra, where low column density DLAs can be identified with more certainty, will be necessary to robustly measure ℓ(X) at z > 4.

Figure 9.

Open in new tabDownload slide

The DLA incidence rate ℓ(X). This observable is more sensitive to the lowest |${\it N}_{\rm H\,{\small i}}$| DLAs than |$\Omega_{\rm H\,{\small i}}$|⁠. Grey points show the uncorrected GGG measurement, and black squares with the corrections applied. Each panel shows a different correction, using either mock or high-resolution spectra for two different authors. There are systematic differences comparable to the statistical errors, and they depend on whether the mocks or the high-resolution spectra are used to estimate correction factors. Similar differences are also found between the two different authors who searched for DLAs. These illustrate that significant systematic uncertainties affect the measurement of ℓ(X).

We can still make a more robust measurement of |$\Omega_{\rm H\,{\small i}}$|⁠, however, regardless of the uncertainty in ℓ(X), as Fig. 10 illustrates. DLAs with the largest contribution to |$\Omega_{\rm H\,{\small i}}$| have |${\it N}_{\rm H\,{\small i}}$| in the range 10^20.8–10^21.6 cm⁻², and DLAs with lower |${\it N}_{\rm H\,{\small i}}$| make a substantially smaller contribution. Therefore, while systematic effects may give rise to a large uncertainty in the number of low column density systems (and thus ℓ(X)), |$\Omega_{\rm H\,{\small i}}$| can still be measured accurately. This point is discussed further in Section 5.

Figure 10.

Open in new tabDownload slide

The differential |$\Omega_{\rm H\,{\small i}}$| distribution, |${\rm d}{\rm \Omega _{\rm H\,\small {I}}}/({\rm d}X{\rm d}\log {\it N}_{\rm H\,{\small i}})$|⁠. DLAs with the largest contribution to |$\Omega_{\rm H\,{\small i}}$| have |${\it N}_{\rm H\,{\small i}}$| in the range 10^20.8–10^21.6 cm⁻²; DLAs with lower |${\it N}_{\rm H\,{\small i}}$| are less important. Therefore, while there is uncertainty in the number of low column density systems (and thus ℓ(X)), |$\Omega_{\rm H\,{\small i}}$| can still be measured accurately. This is also illustrated by Fig. 11.

5 RESULTS AND DISCUSSION

5.1 |$\Omega_{\rm H\,{\small i}}$| measurement

We can now use the |${\it N}_{\rm H\,{\small i}}$|-dependent correction factor k estimated in the previous section to find f_DLA and thus |$\Omega_{\rm H\,{\small i}}$|⁠. For the GGG sample, we count the number of DLAs in a given absorption path, giving each DLA a weight |$k({\it N}_{\rm H\,{\small i}})$|⁠, where |$k({\it N}_{\rm H\,{\small i}})=k_\mathrm{found}({\it N}_{\rm H\,{\small i}})/k_\mathrm{real}({\it N}_{\rm H\,{\small i}})$|⁠. |$k({\it N}_{\rm H\,{\small i}})$| is then estimated as the ratio of the |$\log {\it N}_{\rm H\,{\small i}}$| histograms shown in Figs 5 and 6, with the uncertainty on each bin given by the uncertainties in k_found and k_real added in quadrature.

There are two main contributions to the final error on |$\Omega_{\rm H\,{\small i}}$|⁠. The dominant contribution is the statistical error due to the finite sampling of DLAs: there are 25–30 DLA candidates in each redshift bin, dependent on whether NHMC or JXP's results are used. We estimate this error using 1000 bootstrap samples from the DLA sample. The second is the systematic uncertainty in the correction factor, |$k({\it N}_{\rm H\,{\small i}})$|⁠. We estimate the effect of this uncertainty using a Monte Carlo technique. |$\Omega_{\rm H\,{\small i}}$| is calculated 1000 times, each time drawing |$k({\it N}_{\rm H\,{\small i}})$| from a normal distribution with a mean given by the |$k({\it N}_{\rm H\,{\small i}})$| histogram bin value and σ determined by the uncertainty on that bin, assuming no correlation between uncertainties in adjacent bins. Then the final error in |$\Omega_{\rm H\,{\small i}}$| is given by adding these two uncertainties in quadrature. We confirmed that |${\it N}_{\rm H\,{\small i}}$| error of each DLA (0.2 dex, see Section 4.3) has a negligible contribution compared to these statistical and systematic uncertainties. We also check that using |${\it N}_{\rm H\,{\small i}}$| measurements from the high-resolution spectra, where available, does not significantly change |$\Omega_{\rm H\,{\small i}}$|⁠.

Since we have separate estimates of |$k({\it N}_{\rm H\,{\small i}})$| from the mocks and high-resolution sample, and two authors performed these estimates, we can make four different measurements of |$\Omega_{\rm H\,{\small i}}$|⁠. We use these to gauge the effect on |$\Omega_{\rm H\,{\small i}}$| of estimating corrections from the mocks versus the high-resolution sample, or of any differences in the way the two authors identified DLAs. The results are shown in Fig. 11. The differences between the mocks compared to the high-resolution sample, and between the two authors, are significantly smaller than the uncertainty on any individual |$\Omega_{\rm H\,{\small i}}$| measurement. Therefore, we conclude that neither the methods we use to estimate |$k({\it N}_{\rm H\,{\small i}})$| nor any differences in DLA detection between methods contribute a significant uncertainty to the final |$\Omega_{\rm H\,{\small i}}$|⁠. We caution that this conclusion only holds for the sample of spectra we analyse. New tests of systematic effects may be required for measurements of |$\Omega_{\rm H\,{\small i}}$| using larger samples of DLAs, or using different resolution or S/N QSO spectra.

Figure 11.

Open in new tabDownload slide

|$\Omega_{\rm H\,{\small i}}$| measured by the two authors using the high-resolution sample (top) and mocks (bottom). The uncertainties on |$\Omega_{\rm H\,{\small i}}$| introduced by any differences in selecting DLA candidates between the authors, or between using mocks or the high-resolution sample, are much smaller than the errors shown, which are a combination of the statistical error and uncertainty in the correction factor (see Section 5).

For the remainder of the paper, we use the measurement of |$\Omega_{\rm H\,{\small i}}$| derived using k from the higher resolution sample and measured by author JXP, which is shown in the top-right panel of Fig. 11. This measurement and the 68 per cent confidence interval are given in Table 4. We assume a 20 per cent contribution to |$\Omega_{\rm H\,{\small i}}$| from systems below the DLA threshold, as described in Section3.2.

5.1.1 Is there a bias from gravitational lensing?

5.2 Comparison with previous measurements

Several groups have made measurements of |$\Omega_{\rm H\,{\small i}}$| at z > 4.5 using DLA surveys (Péroux et al. 2003; Guimarães et al. 2009; S10). These are cumulative results – |$\Omega_{\rm H\,{\small i}}$| measurements from each new QSO sample are combined with older |$\Omega_{\rm H\,{\small i}}$| measurements which used a different DLA survey. While combining results in this way maximizes the statistical S/N of the final result, it results in a heterogeneous sample of quasar spectra with different data quality and different DLA identification methods. As shown in Sections 4 and 5.1, at z > 4.4 different identification methods can produce a systematic uncertainty in |$\Omega_{\rm H\,{\small i}}$| which, although smaller than the statistical uncertainties for our current DLA sample, may still be considerable. Since these analyses did not use mock spectra to explore systematic effects, it is difficult to estimate the true uncertainty in |$\Omega_{\rm H\,{\small i}}$| when combining heterogeneous quasar samples with different selection criteria. In contrast, our sample has homogeneous data quality, QSO selection method and DLA identification procedure, and we use mock spectra to test any systematic effects.⁴

Fig. 12 shows our new results together with previous measurements of |$\Omega_{\rm H\,{\small i}}$|⁠, converted to our adopted cosmology. When multiple measurements of |$\Omega_{\rm H\,{\small i}}$| have been made using overlapping QSO samples and the most recent measurement uses a superset of previous QSO samples, only the most recent measurement is shown. For example, the results of S10 include most of the quasars used by Péroux et al. (2003) and Guimarães et al. (2009), so we show only the S10 result. In all such cases, the most recent measurement is consistent with earlier results. Where previous DLA surveys have quoted |$\rm \Omega ^{\rm DLA}_{\rm H\,\small {I}}$|⁠, we convert to |$\Omega_{\rm H\,{\small i}}$| using the relationship |$\rm \Omega _{\rm H\,\small {I}}= 1.2 \rm \Omega ^{\rm DLA}_{\rm H\,\small {I}}/ 1.3$|⁠. Our measurement at 〈z〉 = 4 is higher than, but consistent with, earlier measurements by S10. As such and because we find a possible systematic increase in |$\Omega_{\rm H\,{\small i}}$| towards bright QSOs, we checked whether the magnitude distribution of the S10 QSOs was lower than the GGG sample. z-band data were not available for the whole S10 sample, but the eight QSOs which overlap between their sample and ours have a similar fraction of QSOs with z ≤ 19.2 and z > 19.2 mag. Therefore, a difference in QSO magnitudes is unlikely to cause a difference between our result and the S10 result, and it seems more likely that the difference is caused by a statistical fluctuation.

Figure 12.

Open in new tabDownload slide

Measurements of |$\Omega_{\rm H\,{\small i}}$| at different redshift, from Zwaan et al. (2005), Rao, Turnshek & Nestor (2006), Lah et al. (2007), Braun (2012), Martin et al. (2010), Noterdaeme et al. (2012), Rhee et al. (2013), Delhaize et al. (2013) and S10 (see Table 5). We do not show the measurement using SDSS QSOs by Prochaska & Wolfe (2009); it is consistent with the measurement by Noterdaeme et al. (2012), who use a superset of SDSS QSOs. We also do not show the Péroux et al. (2003) and Guimarães et al. (2009) results, which have a large overlap with the QSO sample used by S10 and are consistent with that measurement. Finally, for clarity, we do not show the measurements at lower redshift from Freudling et al. (2011) and Meiring et al. (2011); they are consistent with the plotted values. All measurements have been converted to the same cosmology (h = 0.7, Ω_m = 0.3, Ω_Λ = 0.7) and include H i mass only, with no contribution from helium or molecular hydrogen.

Our results at 〈z〉 = 4.9 give the most robust indication to date that there is no strong evolution in |$\Omega_{\rm H\,{\small i}}$| over the ∼1 Gyr period from z = 5 to 3. We see a slight drop in |$\Omega_{\rm H\,{\small i}}$| between our z ∼ 4 and z ∼ 4.9 |$\Omega_{\rm H\,{\small i}}$| measurements, but this difference is not statistically significant. If the metal content of DLAs does change suddenly at z = 4.7, as suggested by Rafelski et al. (2014), there is no evidence that it is accompanied by a concomitant change in |$\Omega_{\rm H\,{\small i}}$|⁠. However, the uncertainties remain large and future observations should continue to test this possibility.

Fig. 12 also shows a power law with the form |$\Omega _{\rm H\,\small {I}} = A(1+z)^\gamma$| fitted to the binned data. This simple function provides a reasonable fit (χ² per degree of freedom =1.44) across the full redshift range, with best-fitting parameters A = (4.00 ± 0.24) × 10⁻⁴ and γ = 0.60 ± 0.05. There is no obvious physical motivation for this relation nor any expectation that it should apply at redshifts >5. Nevertheless, it may provide a useful fiducial model to compare to simulations and future observations.

We also compare our new high-redshift value to lower redshift |$\Omega_{\rm H\,{\small i}}$| measurements. As previous authors have noted (e.g. Prochaska et al. 2005; Noterdaeme et al. 2009; Prochaska & Wolfe 2009), |$\Omega_{\rm H\,{\small i}}$| evolves from z = 3 to 0 by factor of ≲ 2, at odds with the very strong evolution in the star formation rate over the same period. Moreover, the drop in |$\Omega_{\rm H\,{\small i}}$| is much smaller than the increase in stellar mass over this period. Fig. 13 demonstrates this point by showing the increase in comoving mass density in stars from z = 5, ρ_⋆ − ρ_⋆(z = 5) and the contemporaneous decrease in H i comoving gas mass density,⁵|$\rho ^{\rm H\,{\small i}}_g(z=5)-\rho ^{\rm H\,{\small i}}_g$|⁠, using the power-law fit from Fig. 12. The mass in stars is calculated using the expression from Madau & Dickinson (2014), and the range shows an uncertainty of 50 per cent, indicative of the scatter in observations around this curve. While the evolution of |$\Omega_{\rm H\,{\small i}}$| from z = 5 to 3 remains uncertain, the H i phase at z = 5 contains ample mass density to form all the stars observed at z ∼ 3, and the evolution predicted by the simple power-law function is consistent with this scenario. From z ∼ 3 to z ∼ 0, however, there is a factor of 5–6 shortfall in H i mass density compared to amount needed to produce stars over the same period. This underscores that at z ≲ 3, the H i phase must be continually replenished by more highly ionized gas, presumably through a combination of cold-mode accretion (e.g. Dekel et al. 2009) and recycled winds (e.g. Oppenheimer et al. 2010). The more highly ionized Lyman limit systems and sub-DLAs should then be important tracers of the interface between this H i phase and more highly ionized gas (e.g. Fumagalli et al. 2011).

Figure 13.

Open in new tabDownload slide

The increase in comoving stellar mass density from z = 5 to 0 (from Madau & Dickinson 2014, thin line and shading) and the corresponding decrease in H i gas mass density over the same period (thick line) using the fitting formula from Section 5.2. Before z ∼ 3, the H i gas phase contains ample mass density to fuel all the observed star formation. However, from z ∼ 3 to the present, it contributes less than ∼20 per cent of the mass necessary to form stars, and so must be continually replenished by more highly ionized gas.

There are several reasons to expect the neutral fraction of the Universe to evolve at z > 3. As we approach the epoch of reionization, the filling factor of neutral hydrogen in the Universe should increase, as large pockets of the Universe are no longer ionized. This is reflected in the decrease in the mean free path for H-ionizing photons (Fumagalli et al. 2013; Worseck et al. 2014) towards higher redshifts. While the bulk of reionization is thought to occur at z > 6, large neutral regions may persist to lower redshifts (e.g. Becker et al. 2015). Our results suggest that while regions of this kind may exist, they do not change the total neutral gas mass density appreciably from that observed at z ∼ 3. This is consistent with the conclusions of Becker et al., who find that by z = 5 the bulk of IGM absorption is due to density fluctuations instead of large, neutral regions yet to be reionized.

This is perhaps not surprising. The distribution of these neutral pockets depends on the nature of reionization, which may progress from low-density regions to high-density regions (‘outside-in’) or the reverse (‘inside-out’), or some combination of the two (e.g. Finlator et al. 2012). However, favoured scenarios see the highest density regions with Δ ≡ ρ/〈ρ〉 ≫ 100 reionized first, as they are populated by galaxies, believed to be the dominant source of ionizing photons. In this case, neutral pockets will persist only in underdense regions such as filaments or voids, with Δ < 100. At z ∼ 2.5 clustering measurements suggest that most DLAs are found inside haloes with masses 10¹⁰–10¹² M_⊙ (Cooke et al. 2006; Font-Ribera et al. 2012), which have a mean Δ > 100. Therefore, even if large neutral regions do persist to z = 5, they may not occur at cosmic densities high enough to produce strong DLA absorption. The remnants of such regions may be observable as Lyman limit systems however, and so one might expect an increase in their incidence rate towards z ∼ 5, which observations already hint may be the case (Prochaska et al. 2010; Fumagalli et al. 2013). The GGG sample can also be used to measure the LLS incidence rate at z > 4, which we will present in a future work.

5.3 Comparison with theory

In Fig. 14 we show |$\Omega_{\rm H\,{\small i}}$| in comparison to some recent theoretical predictions for its evolution. These are by Lagos et al. (2014) using the semi-analytic galform model (see also Popping, Somerville & Trager 2014; Berry et al. 2014), by Bird et al. (2014) from a simulation using the moving-mesh code arepo, and by Tescari et al. (2009, see also Duffy et al. 2012; Rahmati et al. 2015) and Davé et al. (2013) using smoothed particle hydrodynamics (SPH) simulations. While these models broadly match the slow evolution of |$\Omega_{\rm H\,{\small i}}$| since z ∼ 4, most struggle to reproduce the trend of decreasing |$\Omega_{\rm H\,{\small i}}$| with time (with Tescari et al. being a notable exception). Lagos et al. suggest that their model's underestimation of |$\Omega_{\rm H\,{\small i}}$| at high redshift may be due to more neutral gas being found outside galaxy discs in the early Universe. If this interpretation is correct, then our observations suggest that more than half the neutral gas mass (and more than half of DLAs) are found outside galaxies at z ∼ 5. Alternatively, Davé et al. (2013) show that agreement between their simulations and observations can be improved by assuming that a population of low-mass galaxies, unresolved by current SPH simulations, make a significant contribution to the DLA absorption cross-section at high redshift.

Figure 14.

Open in new tabDownload slide

Measurements of |$\Omega_{\rm H\,{\small i}}$| compared to recent theoretical predictions. For clarity, the mean of measurements at z < 0.2 (the error bar shows the standard deviation) is shown. Lines show predictions from a recent semi-analytic model (Lagos et al. 2014), along with SPH (Tescari et al. 2009; Davé et al. 2013) and moving-mesh (Bird et al. 2014) simulations. All the models have been converted to our adopted cosmology. While they do reproduce the roughly flat evolution of |$\Omega_{\rm H\,{\small i}}$| from z = 5 to 0 (in comparison to the cosmic star formation rate density), they do not match the data across the full redshift range.

It is evident that further improvements are needed to theoretical models to reproduce the evolution of |$\Omega_{\rm H\,{\small i}}$| across the full redshift range. If much of the neutral gas is found in galactic outflows or recycled winds, the sub-grid prescription for outflows in SPH simulations may have a strong influence on the predicted |$\Omega_{\rm H\,{\small i}}$| (e.g. Bird et al. 2014). Furthermore, given the small sizes of DLAs (∼5 kpc; Cooke et al. 2010), it may also be important to correct for any smoothing over small-scale density peaks where DLAs are produced, and account for hydrodynamic instabilities which are not resolved by current cosmological simulations (e.g. Crighton et al. 2015).

6 SUMMARY

We have measured |$\Omega_{\rm H\,{\small i}}$| at 3.5 < z < 5.3 using the GGG survey, a homogeneous sample of 163 QSO spectra with emission redshifts >4.4. All the QSOs were colour-selected from the SDSS and so have a well-understood selection function which is independent of any strong absorption in the QSO spectra. Using a combination of higher resolution spectra of DLA candidates and mock spectra, we explore systematic uncertainties in identifying DLAs due to strong IGM absorption at high redshift and the low spectral resolution of the GMOS spectra. The main conclusions from our analysis are as follows.

• We derive the most precise measurement of |$\Omega_{\rm H\,{\small i}}$| at 〈z〉 = 4.9 to date, with a redshift path length at z > 4.5 a factor of 8 larger than previous analyses. |$\Omega_{\rm H\,{\small i}}$| at z = 4.5 is consistent with the value measured at z = 3–3.5, and there is no evidence that |$\Omega_{\rm H\,{\small i}}$| evolves strongly over the Gyr period from redshift 5 to 3. There is also no evidence for an abrupt change in |$\Omega_{\rm H\,{\small i}}$| between z = 4 and 5, which may be associated with a sudden change in metallicity reported at a similar redshift (Rafelski et al. 2014). However, such a change is not strictly ruled out by the data.

• We quantify and correct for the fraction of spurious DLA candidates, and for any DLAs missed in the low-resolution spectra, using higher resolution and mock spectra. We also estimate the uncertainty in the DLA column densities. For this DLA sample, the uncertainty introduced by these systematic effects on the |$\Omega_{\rm H\,{\small i}}$| measurement is smaller than the statistical uncertainties.

• Using the higher resolution spectra and mocks, we show that the typical uncertainty on the DLA |${\it N}_{\rm H\,{\small i}}$| and redshift is 0.2 dex and 200 km s^{− 1}, respectively. Despite the increased IGM absorption at higher redshifts and the low spectral resolution, we find no strong systematic offset in the estimated |${\it N}_{\rm H\,{\small i}}$| for DLAs either as a function of redshift or |${\it N}_{\rm H\,{\small i}}$|⁠.

• We find an excess in |$\Omega_{\rm H\,{\small i}}$| (30 ± 20 per cent) from the brighter half of our QSO sample compared to the fainter half. This is consistent with similar effects found in previous analyses at z ∼ 2.5, which posited gravitational lensing as a possible explanation. Given the smaller Einstein radius at z = 4.5 compared to z = 2.5, for our sample this effect seems more likely to be caused by a statistical fluctuation. As such it should not significantly bias our result.

• Recent theoretical models do not match the data across their full redshift range (z = 5 to 0). A simple power-law model of the form |$\Omega _{\rm H\,{\small i}} = A(1+z)^\gamma$| with A = (4.00 ± 0.24) × 10⁻⁴ and γ = 0.60 ± 0.05, while not physically motivated, does describe the observations over the entire redshift range.

We thank Marcel Neeleman for comments on an earlier version of this paper, Regina Jorgenson for providing a MagE spectrum of one of the GGG DLAs and the referee for their suggestions. Simeon Bird, Claudia Lagos, Romeel Davé and Edoardo Tescari kindly provided tables of their model predictions. We also thank the late Arthur M. Wolfe for providing unpublished ESI spectra and for early contributions to this work.

Our analysis made use of astropy (Astropy Collaboration 2013), xidl⁶ and matplotlib (Hunter 2007). NC and MM thank the Australian Research Council for Discovery Project grant DP130100568 which supported this work. MF acknowledges support by the Science and Technology Facilities Council (grant number ST/L0075X/1). SL has been supported by FONDECYT grant number 1140838 and received partial support from Center of Excellence in Astrophysics and Associated Technologies (PFB 06).

Based on observations obtained at the Gemini and W. M. Keck Observatories. We wish to acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.