Measuring the HI mass function below the detection threshold

TL;DR

This paper introduces a Bayesian stacking method to measure the HI mass function below detection thresholds, enabling accurate estimation of galaxy HI content at cosmological distances using simulated data relevant for upcoming surveys.

Contribution

The paper presents a novel Bayesian stacking technique that allows direct measurement of the HI mass function below the nominal detection limit, improving constraints on galaxy HI content at higher redshifts.

Findings

The HIMF can be accurately reconstructed down to $M_{HI} = 10^{7.5} M_{\odot}$ in simulations.

Constraints on Schechter function parameters improve with lower noise and larger survey areas.

The method is effective for estimating the HIMF at redshifts up to 0.55.

Abstract

We present a Bayesian Stacking technique to directly measure the HI mass function (HIMF) and its evolution with redshift using galaxies formally below the nominal detection threshold. We generate galaxy samples over several sky areas given an assumed HIMF described by a Schechter function and simulate the HI emission lines with different levels of background noise to test the technique. We use Multinest to constrain the parameters of the HIMF in a broad redshift bin, demonstrating that the HIMF can be accurately reconstructed, using the simulated spectral cube far below the HI mass limit determined by the $5\sigma$ flux-density limit, i.e. down to $M_{\rm HI} = 10^{7.5}$ M$_{\odot}$ over the redshift range $0 < z < 0.55$ for this particular simulation, with a noise level similar to that expected for the MIGHTEE survey. We also find that the constraints on the parameters of the Schechter function, $\phi_{\star}$, $M_\star$ and $\alpha$ can be reliably fit, becoming tighter as the background noise decreases as expected, although the constraints on the redshift evolution are not significantly affected. All the parameters become better constrained as the survey area increases. In summary, we provide an optimal method for estimating the HI mass at cosmological distances that allows us to constrain the HI mass function below the detection threshold in forthcoming HI surveys. This study is a first step towards the measurement of the HIMF at high ($z>0.1$) redshifts.

Tables

Table 1. Table 1 : The input parameters of the HIMF model used for our simulation of an H  i galaxy sample, along with the prior range used when running MultiNest . ϕ ⋆ subscript italic-ϕ ⋆ \phi_{\star} has units Mpc -3 dex -1 , while the M ⋆ , M min , M max subscript 𝑀 ⋆ subscript 𝑀 min subscript 𝑀 max M_{\star},M_{\rm min},M_{\rm max} are in units of M ⊙ . For the input value of the cosmological HI density we integrate the input HIMF over the range 7.5 < log 10 ( M HI / 7.5<\log_{10}(M_{\rm HI}/ M ) ⊙ < 12.5 {}_{\odot})<12.5 and divide by the critical density (see equation ( 4 )).

Parameter	Input	Prior Probability Distribution
${\log }_{10}\left({\varphi }_{\star }\right)$	-2.318	uniform $\in [-5,0]$
${\log }_{10}\left(M_{\star }\right)$	9.96	uniform $\in [7,12]$
$\alpha$	-1.33	uniform $\in [-5,0]$
${\log }_{10}\left(M_{\min }\right)$	7.5	uniform $\in [5,10]$
${\log }_{10}\left(M_{\max }\right)$	12.5	uniform $\in [8,13]$
$\beta$	1.0	uniform $\in [-1,4]$
	4.85
	5.24
${\Omega }_{\mathrm{HI}}\times {10}^4$	5.67
	6.12
	6.56
	6.89

Table 2. Table 2 : The parameters of the busy function for simulating the H  i lines. Note that only parameters b 1 , b 2 , c , x p subscript 𝑏 1 subscript 𝑏 2 𝑐 subscript 𝑥 p b_{1},b_{2},c,x_{\rm p} and w 𝑤 w are randomly sampled, the centroid of the error function, x e subscript 𝑥 e x_{\rm e} , is fixed at the redshift of the source. x p subscript 𝑥 p x_{\rm p} and x e subscript 𝑥 e x_{\rm e} are in units of channel, while the line width has units km s -1 .

Parameter	Meaning	Probability Distribution
$b_1,b_2$	steepness of line flanks	uniform $\in [0.3,0.7]$
$c\times {10}^5$	amplitude of the trough	uniform $\in [0.3,0.7]$
$\left|x_p-x_e\right|$	difference of centroids	uniform $\in [0,10]$
$2\times w$	line width	uniform $\in [55,330]$

Table 3. Table 3 : Reconstructed parameters of the HIMF model from our baseline survey with changing background noise, where we fix log 10 ⁡ ( ϕ ⋆ ) = − 2.318 subscript 10 subscript italic-ϕ ⋆ 2.318 \log_{10}(\phi_{\star})=-2.318 for individual redshift bins.

Parameter	$\times {\sigma }_{\mathrm{ch}}$
${\log }_{10}\left({\varphi }_{\star }\right)$	0.5	-2.318	-2.318	-2.318	-2.318	-2.318	-2.318	-2.338${}_{-0.027}{}^{}{}_{}{}^{+0.026}$
${\log }_{10}\left({\varphi }_{\star }\right)$	1.0	-2.318	-2.318	-2.318	-2.318	-2.318	-2.318	-2.397${}_{-0.042}{}^{}{}_{}{}^{+0.039}$
${\log }_{10}\left({\varphi }_{\star }\right)$	5.0	-2.318	-2.318	-2.318	-2.318	-2.318	-2.318	-2.378${}_{-0.152}{}^{}{}_{}{}^{+0.131}$

Equations

$$M_{\mathrm{HI}}=2.356\times 10^5{D}_L^2(1+z{)}^{-1}S,$$

$$\varphi (M_{\mathrm{HI}},z)=\ln (10){\varphi }_{\star }{\left(\frac{M_{\mathrm{HI}}}{M_{\star }}\right)}^{\alpha +1}{e}^{-\frac{M_{\mathrm{HI}}}{M_{\star }}}(1+z{)}^{\beta },$$

$$\Phi (M_{\mathrm{HI}},z_1,z_2)=\frac{{\int }_{z_1}^{z_2}\frac{\mathrm{d}V}{\mathrm{d}z}\varphi (M_{\mathrm{HI}},z)\mathrm{d}z}{{\int }_{z_1}^{z_2}\frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}z},$$

$${\Omega }_{\mathrm{HI}}=\frac{1}{{\rho }_{\mathrm{c}}}{\int }_{{\log }_{10}(M_{\min })}^{{\log }_{10}(M_{\max })}{M}_{\mathrm{HI}}\Phi (M_{\mathrm{HI}},z_1,z_2)\mathrm{d}{\log }_{10}(M_{\mathrm{HI}}),$$

$$\mathcal{L}\propto \prod _{i=bin}{P}_i(k_i|\mathrm{d}N/\mathrm{d}S),$$

$$P_i(k_i|\mathrm{d}N/\mathrm{d}S)=\frac{{\lambda }_i^{k_i}{e}^{-{\lambda }_i}}{k_i!},$$

$${\lambda }_i$$

$$={\int }_{S_{min}}^{S_{max}}\mathrm{d}S\frac{\mathrm{d}N}{\mathrm{d}S}\frac{1}{2}\left[\mathrm{erf}\left(\frac{S_{m_i}-S}{\sqrt{2}{\sigma }_{\mathrm{n}}}\right)-\mathrm{erf}\left(\frac{(S_{m_i}+\Delta S_{m_i}-S}{\sqrt{2}{\sigma }_{\mathrm{n}}}\right)\right],$$

$$\frac{\mathrm{d}N}{\mathrm{d}S}$$

$$=\frac{\mathrm{d}N}{\mathrm{d}{\log }_{10}(M_{\mathrm{HI}})}\frac{\mathrm{d}{\log }_{10}(M_{\mathrm{HI}})}{\mathrm{d}{M}_{\mathrm{HI}}}\frac{\mathrm{d}{M}_{\mathrm{HI}}}{\mathrm{d}S}$$

$$={\int }_{z_1}^{z_2}\frac{\mathrm{d}V}{\mathrm{d}z}\varphi (M_{\mathrm{HI}},z)\frac{{\log }_{10}(e)}{M_{\mathrm{HI}}}2.356\times 10^5{D}_L^2(1+z{)}^{-1}\mathrm{d}z,$$

Full text

††pubyear: 2019††pagerange: Measuring the H i mass function below the detection threshold–4

1 Introduction

Cold gas in the form of neutral hydrogen atoms (H i) provides the reservoir from which molecules are formed and which subsequently go on to fuel star formation in galaxies, from the distant to the nearby Universe. To understand the whole picture of galaxy formation and evolution, it is crucial to know how the H i gas evolves with cosmic time and as a function of environment.

The preferred approach for tracing the H i gas in the local Universe ($z\approx 0$) is via the direct detection of the neutral hydrogen 21 cm hyperfine emission line. Two prime examples of this are the HI Parkes All-Sky (HIPASS) Survey (Barnes et al., 2001; Zwaan et al., 2005), conducted with the 64m Parkes radio telescope in Australia and the Arecibo Legacy Fast ALFA (ALFALFA) survey (Giovanelli et al., 2005; Martin et al., 2010; Jones et al., 2018) conducted by the Arecibo radio telescope with better sensitivity, angular and spectral resolution compared to the Parkes. However, both of these are limited in terms of redshift and sensitivity and have concentrated on surveying large swathes of the local Universe. Measuring the evolution of the H i content of galaxies with redshift requires receivers that work to lower frequencies, coupled with high sensitivity and preferably with relatively high spatial and spectral resolution to avoid confusion (e.g. Delhaize et al., 2013; Jones et al., 2016).

However, we are entering an era of high sensitivity H i surveys with the completion of a variety of new telescopes, in particular the wide-area surveys opened up with Phased Array Feed technology, principally the Australian Square Kilometre Array Pathfinder (ASKAP; Johnston et al., 2007; DeBoer et al., 2009) and the Netherlands Aperture Tile In Focus (AperTIF; Verheijen et al., 2009; Oosterloo et al., 2009), along with the extremely sensitive Meer Karoo Array Telescope (MeerKAT; Jonas, 2009; Jonas & MeerKAT Team, 2016), which is conducting deeper but narrower blind surveys for H i: The MeerKAT International GHz Tiered Extragalactic Exploration (MIGHTEE) Survey (Jarvis et al., 2016) and the Looking At the Distant Universe with the MeerKAT Array (LADUMA) survey (Holwerda et al., 2012; Blyth et al., 2016). One of the most important goals of these surveys is to help us understand the cosmic evolution of H i in the Universe (e.g. Maddox et al., 2016; Meyer et al., 2015).

However, even with the high sensitivity of MeerKAT, the 21 cm emission signal is so weak that only the most HI-massive galaxies will be directly detected beyond the local (z ¿ 0.1) universe, unless very long integration times are used (e.g. Fernández et al. 2016). A proposed method to go further than the planned survey limits is to rely on stacking techniques, which use the known positions of galaxies from surveys at other wavelengths and allows the measurement of simple quantities, such as the average H i mass (e.g. Rhee et al. 2018 and references therein).

Stacking in the context of H i studies usually involves co-adding the spectral line data at the known locations of many galaxies to improve the signal-to-noise at the expense of information on individual galaxies contributing to the stack (e.g. Delhaize et al., 2013; Rhee et al., 2018). Here, we instead present a Bayesian technique to model the distribution of H i masses and not just co-add or average the data. This technique is based on a maximum-likelihood approach (Mitchell-Wynne et al., 2014), which was extended into a fully Bayesian framework by Zwart et al. (2015), who used it to model the source counts of faint radio continuum sources below the nominal detection threshold. Henceforth we simply describe this as Bayesian stacking, although we note that it is not strictly ”stacking” in the normal sense of the word used in the literature.

Specifically, in this paper we extend this technique to determine the H i-mass function (HIMF) directly from the H i-line flux distribution constructed from extracting integrated fluxes across appropriate numbers of spectral channels at the redshift and position of galaxies selected in another waveband. We note that a similar technique has been used to measure the radio luminosity function of optically selected quasars below the noise (Malefahlo et al., 2019).

The idea behind this technique for H i studies is that assuming we can determine the behaviour of the noise, the individual line-integrated 21-cm fluxes measured at the positions and redshifts of known galaxies can provide more information than just the average H i mass. Accurate redshift information, in addition to positional information, is critical and so several spectroscopic surveys are planned to cover both the LADUMA and MIGHTEE fields (e.g. the Deep Extragalactic VIsible Legacy Survey (DEVILS); Davies et al., 2018). Spectroscopic data from the Galaxy and Mass Assembly (GAMA; Driver et al., 2011) survey will provide acccurate redshift information for the ASKAP Deep Investigation of Neutral Gas Origins (DINGO; Meyer, 2009) survey.

We structure this paper as follows. In Section 2, we describe a source-count model for the H i flux distribution and introduce a Bayesian stacking technique for constraining the HIMF model below the detection threshold. We also describe how we generate the H i galaxy sample given an assumed HIMF model and how this is combined with realistic emission line profiles and input into a noisy spectral cube. In Section 3, we apply our method to the simulated data to test the robustness of our method in measuring the HIMF. In Section 4 we summarise our results, highlight some caveats in relation to simulated versus real data, and outline future work based on the real data from the MIGHTEE survey. We use the standard $\Lambda$CDM cosmology with a Hubble constant $H_{0}=67.7$ km$\cdot$s${}^{-1}\cdot{\rm Mpc}^{-1}$, total matter density $\Omega_{\rm m}=0.308$ and dark energy density $\Omega_{\Lambda}=0.692$ (Planck Collaboration et al., 2016).

2 Method

To simulate an H i cube, we first need to construct the H i flux distribution assuming a source-count model and the expected noise properties for a typical survey.

2.1 HIMF model

Under the usual assumptions (e.g. optically thin gas), the H i mass can be converted to the integrated flux (e.g. Meyer et al., 2017) via

[TABLE]

where the H i mass, $M_{\rm HI}$, is in solar masses, the luminosity distance to the galaxy, $D_{L}$, is in Mpc, and the integrated flux $S$ is in Jy km s*-1*. The ($1+z$) factor is needed when $S$ is expressed in units of Jy km s*-1* rather than Jy Hz.

For the HIMF we adopt a Schechter function model, which has been shown to fit the $z\sim 0$ HIMF from both HIPASS (Zwaan et al., 2005) and ALFALFA (Jones et al., 2018), along with a pure density evolution term to characterise the evolution with redshift. Although a strong evolution of H i mass is not expected from the little information we have, including this term does verify the flexibility of our approach. The specific form we adopt is therefore given by

[TABLE]

where $\phi_{\star}$, $M_{\star}$, $\alpha$, and $\beta$ correspond to the normalization, characteristic mass, faint-end slope and power of the redshift evolution respectively. The HIMF represents the intrinsic number density of galaxies in the Universe as a function of their HI mass at a given redshift. The adopted values of these parameters for our models are given in Table 1.

We then define the volume-weighted average H i mass function over the redshift bin defined by $z_{1}<z<z_{2}$ by

[TABLE]

where ${\rm d}V/{\rm d}z$ is the differential co-moving volume at redshift $z$. The H i mass density of the Universe $\Omega_{\rm HI}$ from $z_{1}$ to $z_{2}$ can then be estimated by integrating the $\Phi(M_{\rm HI},z_{1},z_{2})$, which gives

[TABLE]

where $\rho_{\rm c}$ is the critical density of the Universe at $z=0$.

2.2 Generating the galaxy sample

Here we simulate what we expect from the MIGHTEE survey as an example of a future deep H i survey. Using the HIMF described previously with the parameters shown in Table 1, we generate a mass distribution of simulated galaxies over a volume similar to that expected for the MIGHTEE survey (i.e. $\approx 20$ deg2 over the redshift range $0<z<0.55$). The cube has a spatial resolution of $5$ arcsec and a total of 19,504 26 kHz channels from 913.4 MHz to 1420.5 MHz.

Fig. 1 shows the H i mass against redshift for the generated galaxy sample. As we are aiming to explore the HIMF below the nominal noise threshold, the galaxy sample we use is volume-limited. However, we note that in reality we will require an input catalogue of galaxies in order to know where to extract fluxes from within the data cube, and the input catalogue itself will be flux limited. Therefore, any correlation that exists between the input catalogue selection and the H i mass of the galaxies will effect our ability to measure the HIMF and we return to this in Section 3.

2.3 Simulating H i emission line profiles

Once the galaxy samples are generated, we simulate the flux density, $S_{m}(\varv)$, of H i emission lines using the general form of the ”busy function” with a fourth-degree polynomial trough (i.e. equation (4) in Westmeier et al., 2014). The parameters $b_{1}$, $b_{2}$, $c$, $x_{\rm p}$ and $w$ are randomly sampled to assign a variety of line profiles to the H i emission lines. In Fig. 2 we show the flux density as a function of channel for galaxies of a given mass. The major variations are the steepness of line flanks, amplitude of the troughs, centroid of polynomials and widths of the profile. The difference between centroids of the error function and polynomial indicates the asymmetry of line profile, while the line widths approximate different galaxy inclinations. Although the form of lines changes, the integrated fluxes always keep constant for the sources with the same mass and redshift. The ranges used for these parameters are listed in Table 2. We note that these ranges are set to be fairly broad rather than be realistic in order to prove the robustness of our approach. In particular, a double-horn profile with the line width of 150 km s*-1* for the low HI mass samples may not be expected from observations (e.g. Ianjamasimanana et al., 2012), nevertheless one can find many Gaussian profiles easily from our simulations if looking at the flux density with narrower line widths.

We then generate the Gaussian noise with standard deviation similar to that expected from the MIGHTEE survey. Specifically we adopt $\sigma_{\rm ch}=90$ $\mu$ Jy/channel, which is the expected noise per channel with a channel width of $26$ kHz for the MIGHTEE surveys and $0.5\sigma_{\rm ch}$ as an indication of the improvement possible with a deeper survey, such as LADUMA. We also test the case for a much higher noise of $5\sigma_{\rm ch}$. To determine the total flux density per channel across the spectral range, $S_{m}(\varv)$, we simply add the noise to the profile of the busy function and write the total flux density into the simulated cube. We note that only point sources are simulated in this work since our approach is aimed towards the high redshift regime where unresolved sources are expected to dominate given the spatial resolution of current telescopes.

In Fig. 3, we show the simulated emission line as a function of the channel with $0.5\mathrm{\sigma_{ch}}$ and $5\sigma_{\rm ch}$ Gaussian noise. We do not know a priori the shape of the line profile, we therefore fix the line limits that we measure the flux over at a velocity width $N_{\rm ch}\times d\varv=406$ km s*-1* which corresponds to $N_{\rm ch}=74$ at $z=0$ and $N_{\rm ch}=48$ at $z=0.55$. This ensures that we fully cover the range of expected line profiles. The integrated flux is then given by $S_{m}=\sum_{N_{\rm ch}}S_{m}(\varv)d\varv$, where $d\varv=5.5$ km s*-1* at $z=0$. We note that in our idealised case of uniform Gaussian noise $\sigma_{\rm ch}$, the uncertainty associated with summing over all channels is given by $\sigma_{\rm n}=\sqrt{N_{\rm ch}}\sigma_{\rm ch}d\varv$ and is dependent on redshift (i.e. the $d\varv$ increases as $1+z$ and therefore the number of channels required to sample 406 km s*-1* at the redshift of the source decreases as $1+z$ as shown in Fig. 2, thus the noise increases at $\sqrt{1+z}$). This is reasonable for our simulated cube as we already know that the noise follows a Gaussian distribution, but it might be problematic when we deal with real data, where the noise is likely to vary as a function of spectral window, and in such cases we would have to fit the noise distribution as a function of redshift. Furthermore, we would also not expect the signal-to-noise to be constant across the primary beam and return to this point in Section 4.

With the known positions and redshifts from deep optical and/or near-infrared spectroscopy we are able to extract the flux across the number of channels that we expect to encompass the full H i emission line and determine the source counts, which we subsequently model. Note that in this study we do not consider the effects of confusion since the new generation of deep H i surveys will achieve $5-15$ arcsec spatial resolution and have a spectral resolution of a few 10 km s*-1*. Indeed, if we consider how many ”volumetric beams” (or voxels; where in this case a voxel is the beam area multiplied by the number of channels that we integrate over for the H i line) per source we have given our choice of evolving mass function, then we estimate there are around 200 independent ”voxels” per source, which is well above the nominal confusion level. We also note that H i galaxies are not a strongly clustered population (e.g. Papastergis et al., 2013), and assume that lower-mass H i galaxies are actually less clustered in general. As such our simulation of randomly positioned galaxies is probably not far from reality and would still be affected by confusion if it were an issue. This would result in us not accurately recovering the parameters of the input HIMF. However, we do note that close pairs of galaxies in groups will inevitably lead to some small level of confusion, but given the numbers involved this will be an insignificant perturbation on the derived parameters.

We also note that our method could be extended to incorporate confusion noise, either by building it into the model itself (e.g. Chen et al., 2017), or via alternative methods of measuring the noise properties of the spectral cubes at the position of each source (Section 4).

2.4 The Likelihood Function

Based on Bayes’ theorem, the probability of the model given the data is proportional to the likelihood function (i.e. probability of the data given the model). With the assumption that the number of sources in different flux bins are independent, the likelihood for all the flux bins is given by

[TABLE]

where $P_{i}(k_{i}|dN/dS)$ is the probability of obtaining $k_{i}$ objects in bin $i$ given a source-count model for the H i flux distribution ${\rm d}N/{\rm d}S$ (see e.g. Mitchell-Wynne et al., 2014; Zwart et al., 2015).

The measured flux from a H i emission line, $S_{m}$, is a combination of the intrinsic flux from the source, $S$ plus the contribution from the noise, $S_{n}$, which can obviously be positive or negative.

With the assumption that the number of sources in a given patch of sky follows a Poisson distribution, the probability of finding $k_{i}$ objects over the observed sky area in the $i^{th}$ flux bin [$S_{m_{i}},S_{m_{i}}+\Delta S_{m_{i}}$] is given by

[TABLE]

where ${\rm d}N/{\rm d}S$ is our source-count model for the H i galaxies in the observed area and $\lambda_{i}$ is the theoretically-expected number of H i sources in that bin.

If the noise follows a given distribution $P_{n}(S_{n})$, we can write $\lambda_{i}$ as the convolution of $dN/dS$ and $P_{n}(S_{n})$

[TABLE]

where $S_{\rm min}$ and $S_{\rm max}$ are the minimum and maximum flux that can be measured from the data. These correspond to $M_{\rm min}$ and $M_{\rm max}$ in our model via equation (1).

For the idealised case of where the noise follows a Gaussian distribution centered at zero with a variance $\sigma_{\rm n}^{2}$, equation (7) becomes

[TABLE]

where $\sigma_{\rm n}=\sqrt{N_{\rm ch}}\sigma_{\rm ch}d\varv$, and $N_{\rm ch}$ is the number of channels over the flux density profile of the H i line for each galaxy, $\sigma_{\rm ch}$ is the RMS per channel and $d\varv$ is the velocity width of each channel.

We then relate the source-count model with the redshift dependent HIMF, $\phi(M_{\rm HI},z)$, by

[TABLE]

where $z_{1}$ and $z_{2}$ indicate the redshift limits for a given redshift bin and ${\rm d}V=r^{2}\Omega{\rm d}r$ where $r$ is the co-moving distance and $\Omega$ is the solid angle of a given survey.

2.5 Parameter estimation

Here we describe our approach to constrain the parameters of the input mass function so that we compare the best-fitting parameters with the mass function parameters in order to test if we could recover them.

2.5.1 Priors

Priors are the prior knowledge or limits of a parameter before any relevant evidence is taken into account. In Bayesian statistical inference, they provide the sampling parameter space for computing the posterior probability. We assign a uniform prior probability distribution to the power terms $\alpha$ and $\beta$ and adopt uniform logarithmic priors for $\phi_{\star}$, $M_{\star}$, $M_{\rm min}$ and $M_{\rm max}$. These are listed in Table 1.

2.5.2 Multinest

In order to determine the posterior probability distribution and fit the model, we use Multinest111https://github.com/JohannesBuchner/PyMultiNest to sample the prior parameter space (Feroz et al., 2009; Buchner et al., 2014). Multinest is an efficient and robust Bayesian inference tool based on a nested sampling technique (Skilling, 2004), which produces the posterior samples from distributions with an associated error estimate and the Bayesian evidence term allowing model selection. In this paper we accept the median of posterior samples as the best estimate of each parameter considering the complex nature of the posterior, and the set of parameters with the maximum likelihood as the best fitting HIMF. We do not compare different models, we therefore do not use the evidence term. However, we note that an obvious use of the evidence would be to determine the model which best describes the form of the evolution of the HIMF, and we leave this to future work using real data.

2.6 Survey Parameters

For our baseline survey we assume a 1 deg2 survey with noise of $\sigma_{\rm ch}=90$ $\mu$Jy/channel over the redshift range $0<z<0.55$. We then consider the effects of reducing and increasing the rms of the spectral line cube to investigate how the noise properties affect our ability to measure the HIMF accurately. Furthermore, increasing the background noise with a fixed extraction window is also similar to the situation where only photometric redshifts are available, as one has to extract a larger window for fully covering the line profile and this will inevitably integrate more noise into the measured flux. In real observations we are likely to be limited by an optical flux limit for the spectroscopic sample that is required in order to know where to extract fluxes in the spectral line cube. In reality, the supporting spectroscopy will be limited by a complex optical flux limit that is dependent on galaxy colour, morphology and emission line properties. For our purposes, we assume that this flux limit roughly corresponds to a stellar mass limit (although we note that we do not expect this to be a clean one-to-one relation) and many studies have shown that stellar mass and H i mass are correlated, particularly at $M_{\rm HI}<10^{9}$ M⊙ (e.g. Fig. 1 in Maddox et al., 2015). We therefore also investigate how a higher minimum H i mass for our simulated galaxy sample may affect the recovered HIMF. Finally, we also investigate how the survey area influences the constraints by considering a range of surveys from $0.3-20$ deg2.

3 Results

We start by looking at the reconstructed parameters of our baseline survey and then investigate the effects of increasing noise, minimum mass and survey area. Note that we include all simulated samples (i.e. $M_{\rm HI}>10^{7.5}$ M*⊙*) across the redshift range $0<z<0.55$ in the following analysis for exploring the effects of increasing noise and survey area, except for the effect of a higher minimum mass cut.

For analysing the redshift dependence of our approach, we simply split the samples into six individual redshift bins (z-bins) as shown on the top of Fig. 1. Although there should be an optimal way of setting the bins (for instance, the same way as we do for binning the flux below) our current results will only be improved by a better redshift binning scheme.

To capture the shape of the flux distribution in an individual redshift bin accurately, it is sensible to set the flux bins to be uniform in log-flux space since the low flux sources dominate. However, considering that nearly half of low flux sources have negative measured fluxes (due to the noise), we instead apply the Bayesian Blocks method (Scargle et al., 2013) to determine suitable flux bins. The Bayesian Blocks method requires the minimization of a cost function across the data set and allows varying bin widths while being robust to statistical fluctuations, which is ideal for dynamic distributions of the fluxes. The user can determine how many bins are generated through the correct detection rate $p_{0}$ defined by equation (11) of Scargle et al. (2013). A small $p_{0}$ will be more robust to statistical fluctuations in the data, but could be overly coarse. The optimal value of $p_{0}$ is typically determined empirically. At $z>0.1$ we accept the default value of $p_{0}=0.05$, but at $z<0.1$ we use $p_{0}=0.5$ to avoid coarseness of the binning. In general, our results are insensitive to a large range of reasonable values for $p_{0}$.