An HI intensity mapping survey with a Phased Array Feed

Information on matter distribution in the universe is imprinted on the Large Scale Structure (LSS) traced by galaxies. So, measurements of LSS provide essential constraints on cosmological models. Previously, most LSS studies have been made based on optical redshift and lensing surveys such as Galaxy Redshift Survey (2dFGS, Colless et al. 2001), Sloan Digital Sky Survey (SDSS, York et al. 2000), WiggleZ (Drinkwater et al. 2010), Baryonic Oscillation Spectroscopic Survey (BOSS, Anderson et al. 2012) and Dark Energy Survey (DES, Flaugher 2005). Planned future optical surveys include Large Synoptic Survey Telescope (LSST, Ivezic et al. 2008), Euclid (Cimatti & Scaramella 2012) and Dark Energy Spectroscopic Instrument (DESI, DESI Collaboration et al. 2016). However, alternative techniques are also available through radio continuum and spectralline surveys, particularly with the future Square Kilometre Array (SKA) and its pathfinders (e.g., Prandoni & Seymour 2015; Norris et al. 2011; Staveley-Smith & Oosterloo 2015; Meyer 2009; Duffy et al. 2012). These surveys have the potential to detect all massive galaxies at z < 2, and therefore to measure LSS with minimal cosmic variance.

At the current time, the ability to detect individual galaxies in neutral hydrogen (H i) at cosmological distances is limited because of the weakness of the 21-cm emission line with which the mass of H i in galaxies is measured. The most distant galaxy so far detected in the 21-cm H I emission line is at at z = 0.37 detected in the COSMOS H I Large Extragalactic Survey (CHILES) using the Karl G. Jansky Very Large Array (JVLA) (Fernandez et al. 2016). However, the techniques of 'stacking' and 'intensity mapping' (IM), whereby statistical measurements are made by summing the emission from a large sample of galaxies (Rhee et al. 2016, 2018), or by measuring the correlation function in 21-cm images (Chang et al. 2010; Masui et al. 2013), can in principle allow the detection of LSS at redshifts of about unity, even with current radio telescopes. Telescopes such as BAO from Integrated Neutral Gas Observations (BINGO, Battye et al. 2012), Canadian Hydrogen Intensity Mapping Experiment (CHIME, Newburgh et al. 2014), Hydrogen Intensity and Real-time Analysis eXperiment (HIRAX, Newburgh et al. 2016) and Tianlai (Chen 2012) are currently being constructed with the aim of measuring LSS with high accuracy via H I IM. We refer the reader to Bull et al. (2015) for a comprehensive introduction to IM.

In addition to the study of LSS, a further application for IM is to measure the cosmic neutral hydrogen mass density (Ω_HI) at intermediate redshifts. At 1.5 < z < 5, Ω_HI has been estimated using Damped Lyman-α Absorption (DLA) systems which are defined by deep absorption lines in the spectra of quasi-stellar objects (QSOs), due to intervening clouds of hydrogen (e.g., Noterdaeme et al. 2012; Crighton et al. 2015; Bird et al. 2017). At lower redshifts of z < 1.5, Lyman-α absorption lines are shifted to ultraviolet wavelengths and become undetectable using ground-based telescopes. However, at these intermediate redshifts, by cross-correlation between the H I density field and the optical density field in regions of the sky with suitably deep redshift surveys, IM can be applied to infer the amplitude of fluctuations in the H I density field, and thereby the cosmic H I mass density Ω_HI.

Several pioneering experiments in IM have been undertaken. Pen et al. (2009) reported a convincing crosscorrelation signal between H I Parkes All Sky Survey (HIPASS, Barnes et al. 2001) data and optical data from the 6dF Galaxy Survey (6dFGS, Jones et al. 2004, 2009). Chang et al. (2010) reported the first detection of a signal at the cosmologically significant redshift z ∼ 0.8 utilising the Green Bank Telescope (GBT), which was followed by Masui et al. (2013) and Switzer et al. (2013). Anderson et al. (2018) then explored the effect of environment on H I content by cross-correlating Parkes data with the 2dFGRS strip across the South Galactic Pole. The first detection of the power spectrum via an interferometer was obtained by Vujanovic et al. (2012) with the Australia Telescope Compact Array. Furthermore Santos et al. (2015) and others have pointed out that the upcoming SKA will be powerful in conducting future IM surveys.

However, many challenges remain in IM experiments. The first is the thermal noise from the telescope which, at z ∼ 1, is predicted to be 5∼6 orders of magnitude greater than the predicted IM temperature fluctuations. The second is contamination from Galactic and extragalactic foregrounds which are ∼4 orders of magnitude larger than the H i signal, and significantly polarised. Finally, radio frequency interference (RFI) from terrestrial sources is prevalent at low frequencies, and can easily exceed any other source of noise by many orders of magnitude. Thermal noise is uncorrelated with time, and can be suppressed in the normal manner by integration over time and sampling a wide range of frequencies. Foregrounds mainly consist of the diffuse synchrotron and free-free emission from Galactic and extragalactic sources, which are featureless and can therefore be removed by smoothly varying functions if accurate bandpass calibration is available (Condon 1992). More sophisticated foreground subtraction techniques have also been developed for dealing with IM data and similar data from epoch of reionization (EOR) experiments (Wang et al. 2006; Liu & Tegmark 2011; Chapman et al. 2012; Switzer et al. 2013; Wolz et al. 2015). RFI is normally mitigated by flagging and removal in spectral space, although sophisticated adaptive filtering techniques can also provide additional suppression (Briggs et al. 1997; Reynolds et al. 2017).

In this paper, we introduce a new IM survey aiming to measure the large scale H i emission in the redshift range 0.3 < z < 1. This survey was conducted with Parkes telescope equipped with a novel Phased Array Feed (PAF) (Chippendale et al. 2016), which offers several advantages over normal telescope feeds (Reynolds et al. 2017; Deng et al. 2017). The most important advantages of the PAF for the present work are the ability to form multiple simultaneous beams, to cover a large bandwidth, and to present a well-calibrated bandpass without problems associated with spectral ripple. We cross-correlated the measured H i density field with optical density field derived from the catalogues of WiggleZ Dark Energy Survey and detected the cross-correlation signals at redshift z ∼ 0.75.

The paper is structured as follows: in Section 2 we introduce our observational strategy and data reduction pipeline in detail. In Section 3 we verify the pipeline employing simulations. Section 4 details the foreground removal method that we applied to alleviate residual foreground contamination. The final results and estimated cosmological parameters are reported in Section 5, and we summarise in Section 6. Throughout, the cosmological parameters are given by Komatsu et al. (2009), with Ω_m h² = 0.1358, Ω_b = 0.0456, Ω_Λ = 0.726 and h = 0.705.

2.1. H i Observations

The H i data were obtained from observing program P913 (PI: Staveley-Smith), taken using the CSIRO Parkes 64-m telescope in September and October 2016. The observations were conducted utilising CSIRO-built PAF for the Max-Planck-Institute for Radio Astronomy (MPIfR), which was commissioned at Parkes prior to delivery to the Efflesberg 100-m telescope (see Chippendale et al. 2016 for more details). The data considered here were acquired in the so-called band 1, which covers a frequency range of 700 to 1084 MHz. The surveyed regions were chosen to coincide with the footprints of the 0 h, 1 h, 3 h and 22 h fields of the WiggleZ survey (Drinkwater et al. 2010). During the observation, a beamformer was used to produce 17 discrete beams, each with a measured half-power width of 22.3' at the central frequency we considered in this paper. Due to firmware limitations, only 16 beams could be channelised. The beam spacing pitch was 0.35° and, for the purposes of analysis, the beam shape is assumed to be Gaussian. Figure 1 depicts the footprint of the PAF employed in the observations.

Zoom In Zoom Out Reset image size

The flux density calibrators used for the observations were PKS B1934–638 and Hydra A, depending on the RA of the targeted fields during the observation. To maximise the stability of the data, and to minimise the effect of sidelobe and RFI variation, the science observations were obtained in 'shift-and-drift' mode in which the azimuth and elevation of the telescope were fixed, and the telescope therefore scanned in RA at a fixed declination as the Earth rotated. At the frequency we consider here, RFI is mainly terrestrial in nature and there is more gain to be made by minimising exposure to RFI by drift scanning than can be made in the reduction of 1/f noise by scanning the telescope at rates higher than sidereal. The location of the fields observed at Parkes is shown in Figure 2 and Table 1, and the key observational parameters are listed in Table 2.

Zoom In Zoom Out Reset image size

Table 1. WiggleZ Survey Fields Observed Utilising the Parkes Telescope with the PAF

Field name	Total integration time	Calibrators
3 h	22.6 h	PKS B1934-638, Hydra A
1 h	8.7 h	PKS B1934-638
0 h	15.3 h	PKS B1934-638, Hydra A
22 h	13.3 h	PKS B1934-638

Table 2. Key Observational Parameters

Parameter	Value
Bandwidth	384 MHz
Central frequency	891.5 MHz
Spectral resolution	18.5 kHz
Number of channels	20736
Cycle time	4.5 s
Polarisations	2
Beams	16

2.2. H i Data Reduction

Figure 3 plots the RFI occupancy in the whole frequency range, calculated from the fraction of the data at a given frequency where the system temperature was higher than 600 K during calibrator observations. In this paper, we use the RFI-free frequency range of 800 to 820 MHz corresponding to a central redshift of z = 0.75. Subsequent papers will consider other redshift ranges.

Zoom In Zoom Out Reset image size

The data reduction procedure is as follows:

• We converted the raw data from arbitrary units to units of Kelvin. To do this, we calculated the scaling factors using on-source and off-source spectra for the calibrators, scaling the raw data in the targeted fields with the corresponding frequency-dependent factors. These scaling factors were measured to have uncertainties of ∼ 11%.
• For each scan and each frequency channel, we calculated the running median for ±30 integration cycles, and subtracted the lowest median value from that channel. In the selected frequency band, the system temperatures of some beams were also very high, up to ∼ 800 K, due to bad beam-forming calibration. We discarded beams with a median system temperature higher than 500 K.
• We then subtracted smooth foreground emission and RFI in an iterative way. Given that the foreground is mainly made up of Galactic synchrotron radiation and free-free emission with a monotonic variation with frequency, we removed a linear fit to the spectral data for each integration cycle. Before going to the next step, we iteratively masked 3σ outliers until the result converged.
• To remove long-term temporal variations, we removed a linear fit to the data in each spectral channel for all scans. Then we iteratively masked 3σ outliers until the result converged.
• We calculated how much data had been flagged in each cycle, each channel and each scan. Cycles or channels with high masked fractions compared to the average were completely flagged. Where cycles or channels were bracketed with masked data on either side, the mask was dilated to also mask these cycles or channels. We again iteratively masked 3σ outliers until the result converged.
• The calibration procedures (iii) − (v) were repeated twice.

Figure 4 shows an example of the different steps in the calibration process, with the bottom panel featuring an example of the final calibrated data.

Zoom In Zoom Out Reset image size

The key parameters of the final gridded data cubes are summarised in the first row of Table 3. Grid points were computed using the weighted median of all recorded brightness temperature values (all beams, all times, both linear polarisations) within a radius of half a beam around each given grid point. The weights were computed as the product of two parameters: (1) the reciprocal of the product of root mean square (rms) noise levels in the channel and cycle containing the observed value; and (2) the normalised value of the Gaussian at the position of the observation, where the Gaussian is centred on the grid point and has the same width as the beam. The angular spacing is a quarter of the beam size, 5.58' per pixel, with a frequency spacing of 18.5 kHz, corresponding to a comoving transverse size of 4.3 h⁻¹ Mpc and a comoving line-of-sight distance of 0.11 h⁻¹ Mpc at the central frequency channel.

Table 3. Key Parameters of the Gridded Data Cubes and Simulations in the WiggleZ 3 hr, 1 hr, 0 hr and 22 hr Fields

		Gridded area (J2000)	Volume (Mpc3)	Dimensions
Observation	3 hr field	02h53m < RA < 03h28m, −18°36' < Dec < −05°42'	603 × 412 × 124	(105, 104, 1080)
	1 hr field	00h30m < RA < 01h05m, 00°59' < Dec < 05°12'	200 × 408 × 124	(46, 94, 1080)
	0 hr field	23h21m < RA < 23h58m, −09°00' < Dec < 01°24'	486 × 425 × 124	(112, 98, 1080)
	22 hr field	21h22m < RA < 22h00m, −05°00' < Dec < 04°42'	456 × 451 × 124	(139, 95, 1080)
Simulation	3 hr field	02h08m < RA < 04h13m, −19°21' < Dec < −04°57'	668 × 477 × 124	(120, 119, 1080)
	1 hr field	23h45m < RA < 01h50m, 00°14' < Dec < 05°57'	265 × 473 × 124	(61, 109, 1080)
	0 hr field	22h36m < RA < 00h43m, −09°45' < Dec < 02°09'	556 × 495 × 124	(128, 114, 1080)
	22 hr field	20h37m < RA < 22h45m, −05°45' < Dec < 05°27'	521 × 516 × 124	(154, 110, 1080)

Finally, the rms noise level was calculated for each channel in the gridded data cubes. In each field we iteratively fit a third-order polynomial to the rms spectrum and masked channels with an rms deviating by more than 3σ. As RFI appears in almost the same frequency channels in all fields, we merged the masks from all four fields and applied to the individual fields. Figure 5 depicts the corresponding spectral masks for the individual fields, and the resultant merged mask.

Zoom In Zoom Out Reset image size

The resultant data cube also has different noise levels in different directions. We calculated the rms noise level for each masked spectrum in the data cubes and further masked 22.0% of the highest rms in the areas withWiggleZ coverage (see Appendix). The resultant cubes were downsampled by averaging every 20 channels, resulting in an rms noise level of 35.5 mK for the whole data set. Figure 6 displays the final gridded H i maps we obtained in the four fields at the central frequency channel.

Zoom In Zoom Out Reset image size

2.3. Optical Data

The optical catalogue we relied on for cross-correlation is from the WiggleZ Dark Energy survey which measured redshifts for more than 240 000 galaxies up to redshift z ∼ 1 in seven fields with a total area of approximately 1000 deg². Full details on the design and calibration are given in Drinkwater et al. (2010).

The optical density fields are derived from the correspondingWiggleZ catalogs from the final data release (Drinkwater et al. 2018), which have 1446, 361, 733 and 1407 galaxies within the redshift range of the 3 hr, 1 hr, 0 hr and 22 hr regions, respectively. Due to the optical sample being GALEX-selected, in some areas near bright stars it had no input data.

We gridded the optical data into cubes, representing each galaxy by a 3-dimensional (3D) Gaussian with a full width at half maximum of 150 km s⁻¹ in velocity space and 223 in angular extent (to match the radio data). The units for this spectrum are arbitrary because only the optical overdensity field is required. The exact shape of the spectrum is also arbitrary in that the WiggleZ galaxies are only regarded as tracers of the underlying density field.

To allow for WiggleZ selection effects, we generated 1000 random optical catalogues for which the selection functions and completeness are maintained while positions and redshifts of galaxies are randomised. We computed these randomised optical density fields in the same way and then the real optical density field is corrected as follows

$$\begin{eqnarray}&&{\delta }_{{\rm{opt}}}=D/R\times \langle R\rangle /\langle D\rangle -1\,\,,\end{eqnarray} \tag{ 1 }$$

where D is the original computed optical density field and R is the average density field from random catalogues at the same position. To avoid error introduced by dividing by a small R value, we then masked directions with R < 1, which excludes ∼ 8 per cent of all directions.

Figure 7 shows the original optical density fields, the ⟨R⟩/⟨D⟩ applied for correction and the final corrected and masked optical density fields in four fields.

Zoom In Zoom Out Reset image size

To understand the effect of our pipeline on the observed signal, we simulated an IM survey. In our simulation, we assume that the distributions of neutral hydrogen and dark matter can be represented by the power spectrum generated utilising the CAMB package ¹ (Lewis et al. 2000). To create a simulated data cube, we first created the 3D matter distribution in Fourier space. The values of the voxels in this data cube are complex numbers, with the square of their modulus being related to the value of the power spectrum.

To cover the whole gridded area, we set the size of each simulated field to be 1.5 degrees larger than the real observation along both the RA and Dec directions, and with a spectral resolution the same as the gridded H I cubes. After performing inverse Fourier transform, the comoving volumes of the generated neutral hydrogen density fields are larger than the sizes of the real H I data cubes previously discussed. The lower part of Table 3 summarises the key parameters of the simulations in the four fields.

We convolved the simulated density field with the Gaussian PAF beam pattern and assigned coordinates to each voxel. We then conducted a mock observation of the simulated density field following the actual scanning strategy, assuming identical system temperatures to the real data, and applied our gridding and masking algorithm to the mock data files.

Figure 8 depicts an example of a single plane of the convolved simulated density field and the corresponding mock observation at the central frequency channel. We cross-correlated both, and compared it with the autocorrelation of the former. The result indicates that the observed correlations are 6.3%, 8.1%, 6.7% and 6.8% smaller than the true correlations in the 3 hr, 1 hr, 0 hr and 22 hr field respectively, demonstrating the effect of the observing procedure and pipeline distortions on the signal strength. For the whole masked dataset, an approximate average correction of 7.4% is adopted.

Zoom In Zoom Out Reset image size

Despite careful attention to flagging and calibration, residual foreground and RFI contamination remain in the data at a level that is higher than the H i signal we seek to measure. We therefore applied the technique of singular value decomposition (SVD) to our data cube, one of several foreground removal techniques applied by other authors to IM data sets (Chang et al. 2010; Masui et al. 2013; Anderson et al. 2018). SVD is a generalization of standard eigenvalue decomposition techniques which can be applied to non-square matrices. It decomposes the data cube into orthogonal modes of decreasing amplitude, with the low-order modes mainly containing the unwanted foreground and RFI residuals. In the SVD process, the 3D H i data cubes are first transformed to four 2-dimensional (2D) matrixes, and then are concatenated to an m × n matrix M, where m is the total pixel number on the plane perpendicular to line-of-sight and n is the channel number. The decomposition of matrix M is given by

$$\begin{eqnarray}&&M=U\Sigma {V}^{\star }\,\,,\end{eqnarray} \tag{ 2 }$$

where U is the m × m matrix of left-singular vectors, Σ is an m × n matrix where the singular values are the non-negative reals on the diagonal and V^⋆ is the conjugate transpose of V, which is the n × n matrix of right-singular vectors. The singular values of M represent the amplitudes of each mode after decomposition. As the amplitudes of foreground residuals are usually assumed to be higher than the H i signal, the low-order SVD modes are expected to mainly consist of these residuals, which facilitate a 'cleaner' cross-correlation with the optical density field once these modes are removed and the H i data cube is reconstructed from M. Figure 9 features examples of the recovered zoom-in imaged H i data plane in the 3-hr field at the central channel after the removal of 0, 30 and 120 modes, from top to bottom respectively.

Zoom In Zoom Out Reset image size

The cross-correlation of the H i and optical data is displayed in Figure 10 as a function of the number N of SVD modes removed. The fluctuations in the cross-correlation are large at small N where the foreground and RFI components are dominant. At large N the cross-correlation signal becomes low due to the SVD process of removing genuine H i signal.

Zoom In Zoom Out Reset image size

In order to estimate this signal loss, we combined the optical density field with the H i cube using a scale factor A, replacing ΔT_b with ΔT_b + Aδ_opt in the cross-correlation. As shown in Figure 10, this results in an increased/decreased cross-correlation amplitude which, as expected, also diminishes with N. The loss factor as a function of N can be derived from the ratio of these measurements, as plotted in the lower panel of Figure 10.

The ratio of the measured cross-correlation values to the loss factor therefore gives N estimates of corrected cross-correlation signal. Rather than selecting an arbitrary value for N, we have simply applied a median fit of the loss curve to the measured signal, as represented by the black solid line in the upper panel of Figure 10. This results in a median cross-correlation signal.

To account for signal loss due to the data reduction and gridding process (see Sect. 3), we apply a 7 per cent correction and derive our final result of ⟨ΔT_b δ_opt⟩ = 1.32±0.42 mK. The significance of this result is at the 3σ level, so it is a tentative detection.

The errors represent the standard deviation of the cross-correlation derived from 100 random optical catalogues in which the same optical selection function as the WiggleZ survey is applied, while the positions and redshifts have been randomised. The error is therefore purely statistical.

To estimate errors arising from our processing and analysis pipeline, we tested the pipeline with different parameters, including the width of the velocity convolution function, the upper limit for N used in fitting the correlation amplitudes and the fitting method itself. With reasonable changes, these parameters only cause minor changes (within 10%) to the final results.

However, other more major sources of errors not accounted for include cosmic variance and the discrete sampling of the optical density field. The sparseness of the WiggleZ data in particular may be a major limitation in accurate calculation of the optical density field. The simulations required to model this are outside the scope of this paper.

In principle, the detected signal can be related to the cosmic H i density ΩH i utilising the following relation

$$\begin{eqnarray}&&\langle \Delta {T}_{{\rm{b}}}{\delta }_{{\rm{opt}}}\rangle =44br\langle {\delta }_{{\rm{opt}}}^{2}\rangle \left(\displaystyle \frac{{\Omega }_{{\rm{HI}}}h}{2.45\times {10}^{-4}}\right)\displaystyle \frac{{(1+z)}^{2}}{E(z)}{\rm{\mu}} {\rm{K}}\end{eqnarray} \tag{ 3 }$$

(Chang et al. 2010), where ΔT_b is the H i brightness temperature fluctuation, δ_opt is the normalised optical density field, $b={\langle {\delta }_{{\rm{HI}}}^{2}\rangle }^{0.5}/{\langle {\delta }_{{\rm{opt}}}^{2}\rangle }^{0.5}$ is the bias factor, $r=\langle {\delta }_{{\rm{HI}}}{\delta }_{{\rm{opt}}}\rangle /{(\langle {\delta }_{{\rm{HI}}}^{2}\rangle \langle {\delta }_{{\rm{opt}}}^{2}\rangle )}^{0.5}$ is the stochasticity, E(z) = H(z)/H₀ and h is the Hubble constant in units of 100 km s⁻¹ Mpc⁻¹. However, as the terms in the above are model-dependent, and as there are potentially large uncertainties in correcting for the (sparse) optical density field noted above, realistic galaxy simulations are required, so this derivation is deferred to a later paper.

We have conducted an H i IM experiment with the Parkes 64-m telescope in four equatorial fields which have redshift data out to z ∼ 1 from the WiggleZ redshift survey. The H i data were taken with a novel PAF which permitted simultaneous observations over 16 separate beams. This paper presents the data reduction procedure in a redshift band at z ∼ 0.75 which was judged clear of RFI.

The rms noise level in the data cubes in the range 800–820 MHz was 35.5 mK at a spectral resolution of 0.37 MHz. After foreground removal using an SVD technique, the cross-correlation signal with the WiggleZ optical density fields in the same redshift range was detected at close to the 3σ level. After correction for signal loss caused by the data reduction and analysis pipeline, the amplitude of the observed cross-correlation is ⟨ΔT_b δ_opt⟩ = 1.32 ± 0.42 mK (statistical errors).

Despite the presence of high system temperature and strong RFI, we have been able to make useful IM measurements at a cosmological redshift with a PAF. The interpretation of the cosmic H i density that corresponds to the measured cross-correlation signal (or upper limit) requires a more sophisticated simulation that is beyond the scope of this paper. However, preliminary indications are that a future generation cryogenic PAF will be able to substantially improve the accuracy of these observations and, with more attention to foreground contamination, may even permit a detection of the auto-correlation signal, thus dispensing with the requirement for a matching optical redshift survey.

We thank Chris Blake and Yi-Chao Li for their constructive discussion and suggestions. LL thanks the Chinese Scholarship Council and University of Western Australia for the financial support. We thank the staff of the Parkes radio telescope for their technical assistance. The Parkes radio telescope is part of the Australia Telescope National Facility which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. Parts of this research were conducted by the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020 and the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. This work was also supported by the National Key R&D Programme of China under grant number 2018YFA0404603.

The cross-correlation which we aim to measure can be expressed as

$$\begin{eqnarray}&&C=\langle {\delta }_{{\rm{opt}}}\Delta {T}_{{\rm{HI}},{\rm{o}}}\rangle =\langle {\delta }_{{\rm{opt}}}\Delta {T}_{{\rm{HI}},{\rm{r}}}\rangle +\langle {\delta }_{{\rm{opt}}}{T}_{{\rm{n}}}\rangle,\end{eqnarray} \tag{ A.1 }$$

where ΔT_HI,o, ΔT_HI,r and T_n are the observed H i brightness temperature, real H i brightness temperature and noise, respectively. As the expectation values E(δ_opt) = 0, E(T_n) = 0 and T_n are uncorrelated with δ_opt, the variance of C can be computed as

$$\begin{eqnarray}\begin{array}{lll}D(C) & = & D({\delta }_{{\rm{opt}}}\Delta {T}_{{\rm{HI}},{\rm{r}}})+D({\delta }_{{\rm{opt}}}{T}_{{\rm{n}}})\\ & = & D({\delta }_{{\rm{opt}}}\Delta {T}_{{\rm{HI}},{\rm{r}}})+D({\delta }_{{\rm{opt}}})D({T}_{{\rm{n}}}).\end{array}\end{eqnarray} \tag{ A.2 }$$

The cross-correlation is calculated by averaging from all voxels: ∑C_j /N, where N is the number of voxels used for averaging, and the subscript j denotes the voxel number. The expected variance of the mean estimate of C, assuming unweighted averaging and a universal D(δ_opt), can therefore be approximated as

$$\begin{eqnarray}\begin{array}{lll}\displaystyle \frac{m}{{N}^{2}}\displaystyle \sum _{j=1}^{N}{D}_{j} & = & \displaystyle \frac{m}{{N}^{2}}\displaystyle \sum _{j=1}^{N}{D}_{j}({\delta }_{{\rm{opt}}}\Delta {T}_{{\rm{HI}},{\rm{r}}})\\ & & \,+\displaystyle \frac{m}{{N}^{2}}\displaystyle \sum _{j=1}^{N}{D}_{j}({\delta }_{{\rm{opt}}}){D}_{j}({T}_{{\rm{n}}})\\ & = & \displaystyle \frac{m}{{N}^{2}}\displaystyle \sum _{j=1}^{N}{D}_{j}({\delta }_{{\rm{opt}}}\Delta {T}_{{\rm{HI}},{\rm{r}}})\\ & & \,+D({\delta }_{{\rm{opt}}})\displaystyle \frac{m}{{N}^{2}}\displaystyle \sum _{j=1}^{N}{D}_{j}({T}_{{\rm{n}}}),\end{array}\end{eqnarray} \tag{ A.3 }$$

where m is the number of pixels per area. The first two terms on the right side of this equation are from the real density fields and cannot be reduced. So, the uncertainty may be expected to reach a minimum when the final term S = ∑D_j (T_n)/N² is minimised. In the observation, ΔT_n >> Δ T_HI,r and the final gridded data quality depends much more on the pointing direction rather than frequency. We therefore compute the noise variance for each spectrum, D_j (T_n), and calculate S = ∑D_i (T_n)/N².

In keeping with our preference to mask poor quality data, rather than applying continuous weights, we now eliminate directions with the highest variance. The appropriate metric is to calculate the value for N for which the above quantity S = ∑D_i /N² is minimised after sorting the D_i in ascending order. Figure A.1 shows S a function of N for our data. The result affirms that S reaches a minimum when 22.0% of the directions with the highest rms noise level are masked.

Zoom In Zoom Out Reset image size