2.1. The cosmic gas density

The cosmic gas density, $\Omega_{{\rm H\,I}}(z),$ is a measure of how the global supply of cold gas available for star formation changes with redshift or cosmic time. Whilst the cosmic star formation rate is known to have decreased by an order of magnitude from $z\,{\sim}\,1$ to the present time (Lilly et al. Reference Lilly, Le Fevre, Hammer and Crampton1996; Madau, Pozzetti, & Dickinson Reference Madau, Pozzetti and Dickinson1998), somewhat surprisingly, $\Omega_{{\rm H\,I}}(z)$ seems more constant. Recent semi-analytic (Kim et al. Reference Kim2015) and hydrodynamic (Davé et al. Reference Davé, Rafieferantsoa, Thompson and Hopkins2017) models have now been adjusted to reflect this slower evolution, but there remain modest differences between models at $z<1$ and enormous divergence at $z>1$ (see Figure 19 in Xi et al. Reference Xi2021). Moreover, the observational and systematic errors in measurements of $\Omega_{{\rm H\,I}}(z>0.1)$ remain stubbornly significant. The main methods available to measure cosmic H i density are: (1) direct H i spectral-line measurements of individual galaxies in a blind survey; (2) spectral stacking using optical redshifts; (3) intensity mapping; and (4) damped Lyman- $\alpha$ (DLA), or proxy measurements. Of these, the first is strongly preferred but, owing to the weakness of the 21-cm line, it has been difficult to accomplish at $z>0.1$ . Most of $\Omega_{{\rm H\,I}}$ values at $z > 0.1$ (which are much less accurate than for $z<0.1$ ) use stacking (Lah et al. Reference Lah2007; Delhaize et al. Reference Delhaize, Meyer, Staveley-Smith and Boyle2013; Rhee et al. Reference Rhee2018), intensity mapping (Chang et al. Reference Chang, Pen, Bandura and Peterson2010; Masui et al. Reference Masui2013) or damped Lyman- $\alpha$ (DLA) (Rao & Turnshek Reference Rao and Turnshek2000; Prochaska & Wolfe Reference Prochaska and Wolfe2009; Grasha et al. Reference Grasha, Darling, Leroy and Bolatto2020) methods. With the FUDS, we plan to perform the most accurate measurements to date for $\Omega_{{\rm H\,I}}$ at $z > 0.1$ , without reference to optical surveys. Stacking methods will allow accurate measurements in the range $z=0.2$ – $0.4$ , but with a strong dependence on the selection and assumed evolutionary parameters of the corresponding optical input catalogue.

2.2. The evolution of galaxies and their H i content

The study of the cold gas (and stellar) content of individual galaxies, and the effect of environment on gas supply and star-formation rate, is a key to better understand galaxy evolution. Again, due to the difficulties of individual galaxy detection, this has been challenging. Catinella & Cortese (Reference Catinella and Cortese2015) measured the H i content in 39 massive late-type galaxies at $0.17 \leq z \leq 0.25$ . Compared with scaling relations of galaxies at $z \sim 0$ , they found that these galaxies are the analogues of extreme objects in the local Universe. As above, stacking techniques (Chengalur, Braun, & Wieringa Reference Chengalur, Braun and Wieringa2001; Zwaan, van Dokkum, & Verheijen Reference Zwaan, van Dokkum and Verheijen2001) can shed more light on the mean properties of galaxies over a larger redshift range. For example, Rhee et al. (Reference Rhee2018) measured a higher H i mass-to-luminosity ratio for distant late-type galaxies as compared to samples at lower redshift. However, the requirement to use optical/IR samples adds extra measurement uncertainties. More importantly, it adds selection effects, such as imaging and spectroscopic incompleteness (Jones et al. Reference Jones, Peterson, Colless and Saunders2006), and systematic uncertainties due to redshift-dependent luminosity evolution (McNaught-Roberts et al. Reference McNaught-Roberts2014) and k-corrections (Lake & Wright Reference Lake and Wright2016), and possible environmental dependence (McNaught-Roberts et al. Reference McNaught-Roberts2014). There is also a dependence on optical luminosity and/or luminosity density functions which may have been derived from different optical samples. Perhaps the most important of these additional uncertainties is the tendency for optical catalogues to miss the substantial population of H i-rich but optically low-surface brightness galaxies. This tendency is strongly reflected in the steeper low-mass slope in the H i mass function compared with the optical luminosity or stellar mass functions (Butcher et al. Reference Butcher, Schneider, van Driel and Lehnert2018).

FUDS is intended to be deep enough to detect individual galaxies right up to the survey limit of $z \sim 0.4$ , and in sufficient statistical numbers to allow an accurate study of the evolution of galaxies and their gas content. However, given the large physical size of the FAST beam at higher redshifts, confusion is a possible problem (see e.g. Jones et al. Reference Jones, Haynes, Giovanelli and Papastergis2016a). We will briefly discuss the impact of confusion on our results in Section 7.

2.3. The evolution of the H i mass and velocity functions

As statistical quantities which relate to the number density of galaxies of different gas and halo mass, the HI Mass Function (HIMF) and velocity width function (HIWF) are important constraints for cosmological models. Thanks to the large sample sizes now available, the form of the HIMF and HIWF and their Schechter parameterisation (Schechter Reference Schechter1976) have been greatly refined (Jones et al. Reference Jones, Haynes, Giovanelli and Moorman2018; Oman Reference Oman2022). However, the HIMF beyond the local Universe is poorly known. Using the Arecibo Ultra-Deep Survey (AUDS), Hoppmann et al. (Reference Hoppmann2015) did not find any evolutionary trend in the preliminary AUDS60 sample. However for the full sample, Xi et al. (Reference Xi2021) found weak evidence for a ‘knee’ mass which increases with redshift. This conclusion was drawn based on a sample of only 247 galaxies. With its larger sample size, and wider redshift range, FUDS will put more accurate constraints on evolutionary trends.

2.4. Environment dependence

Many studies have explored the environmental dependence of the HIMF. For example, does the density ratio of low-mass and high-mass galaxies increase or decrease with local density? No common trends have yet been established, either for the slope of the low-mass end of the HIMF ( $\alpha$ ) or the value of the characteristic mass ( $M^*$ ). Rosenberg & Schneider (Reference Rosenberg and Schneider2002) found a flatter $\alpha$ in an intermediate density region around the Virgo cluster compared with the field. Busekool et al. (Reference Busekool2021) also found a flat slope, $\alpha=-0.92$ , in dense regions using H i detections in the Ursa Major cluster. In contrast, Pisano et al. (Reference Pisano2011) found flatter $\alpha$ as the density decrease based on six loose groups. Zwaan et al. (Reference Zwaan, Meyer, Staveley-Smith and Webster2005) found steeper $\alpha$ values in dense regions based on the full HIPASS sample. None of these studies found any change in $M^*$ . Springob et al. (Reference Springob, Haynes and Giovanelli2005) used ${\sim}3\, 000$ AGC galaxies to derive the HIMF and found constant $\alpha$ and decreasing $M^*$ in dense regions. Utilising the 40% and 70% ALFALFA samples, Moorman et al. (Reference Moorman2014) and Jones et al. (Reference Jones, Papastergis, Haynes and Giovanelli2016b) discovered positive correlation between $M^*$ and density, respectively. But, neither found any obvious change in $\alpha$ . Jones et al. (Reference Jones, Haynes, Giovanelli and Moorman2018) derived a steeper $\alpha$ and higher $M^*$ for non-cluster high-density regions from the 100% ALFALFA sample. Said , Kraan-Korteweg, & Staveley-Smith (Reference Said, Kraan-Korteweg and Staveley-Smith2019) found flatter $\alpha$ in dense regions within the Parkes H i Zone of Avoidance survey (HIZOA, Donley et al. Reference Donley2005; Staveley-Smith et al. Reference Staveley-Smith2016). However, they did not find any change in $M^*$ . By examining HI detections, and the form of the HIMF in different target regions, H i surveys such as FUDS, with complementary data to help define density metrics, will provide better insights into the effect of environment.

3.1. FAST

The Five-hundred-metre Aperture Spherical Telescope (FAST, Nan et al. Reference Nan2011) has a spherical main reflector with a diameter of 500 m, which is part of a sphere of radius 300 m. FAST is located at longitude $106^\circ 51^{\prime} 24.\!\!^{\prime\prime} 00074\ \textrm{E}$ , latitude $25 ^{\circ} 39^{\prime} 10.\!\!^{\prime\prime} 62653\ \textrm{N}$ , Observations with zenith angles smaller than $26.4^{\circ}$ can be carried out with full gain (Jiang et al. Reference Jiang2020). In this mode, the accessible sky area is $-0.7^{\circ}$ < Dec (J2000) $< 52.1^{\circ}$ . FAST can also work to a maximum of zenith angle of $40^{\circ}$ with lower gain and higher system temperature. The observable sky area is then enlarged to $-14.3^{\circ} < \textrm{Dec}$ (J2000) $< 65.7 ^{\circ}$ .

During an observation, actuators beneath the main reflector dynamically change the shape of part of the surface to form a parabola with an illuminated diameter of 300 m. The prime focus is used to collect electromagnetic radiation across the focal plane. The design focal ratio was 0.4665, but an optimised ratio of 0.4621 is presently adopted, which gives a focal length of 138.63 m. The position of the dome containing the primary focus and feed system is controlled by six cables, with little blockage of the optical path, resulting in low RFI scattering and low far-sidelobes compared with conventional dishes with massive support struts.

The front-end used for our observations is the 19-beam receiver (Jiang et al. Reference Jiang2020). According to the angular distance from centre beam, the receiver has one central beam (01), 6 beams (02 – 07) evenly spread around an inner circle, 6 beams (09, 11, 13, 15, 17, 19) around a second concentric circle, and 6 beams (08, 10, 12, 14, 16, 18) around the third circle. The inner, second, and third circles are shown by light blue circles in Figure 1. The angular distance between the closest beams is 5.9 arcmin, and the spacing between parallel rows in the geometry shown in Figure 1 is 5.11 arcmin. Each beam has two orthogonal linear polarisation channels (X and Y). The receiver can rotate to track the parallactic angle. The allowable rotation of the line joining beams 08, 02, 01, 05, 14 with a line of constant declination (current epoch) is $\pm 80 ^{\circ}$ .

Figure 1. The geometry of the 19 beams of the FAST multibeam receiver is shown on the lower right corner. The black solid circles show the beam size of 2.9 arcmin at 1.4 GHz. The angular distance between nearby beams is 5.9 arcmin, and the spacing between nearby rows is 5.11 arcmin. The index of each beam is given in the circles. The light blue circles indicate the separation of different beams from the central beam. The black star is the position of the FUDS0 (continuum) calibrator source. The desired trace of this example calibration scan is shown with the black dashed line. The real trace of the telescope is shown with the red dotted line.

A calibration noise diode is installed at $45^{\circ}$ relative to two linear polarisation probes in each horn. The noise diode operates at two temperatures (high at 10 K or low at 1 K) for different requirements. Since the structure of FAST main reflector depends on pointing direction, calibration of main reflector is frequently required.

The FAST backend consists of 12 Re-configurable Open Architecture Computing Hardware version 2 (ROACH2) units developed by the Collaboration for Astronomy Signal Processing and Electronics Research (CASPER; Hickish et al. Reference Hickish2016; Jiang et al. Reference Jiang2019). Each ROACH2 unit records two pairs of polarised signals simultaneously. Hence, 10 are used to record the signals from two polarisations of 19 beams, and 2 are spare in case of breakdown.

The bandwidth of the ROACH2 is configured to be 500 MHz (frequency range 1.0–1.5 GHz), split into a maximum 1 048 576 channels with a maximum frequency resolution of 476.8 Hz. For our observations, we use 65 536 channels, which gives a frequency spacing of 7.63 kHz, equivalent to a velocity spacing of $1.61\ \textrm{km} \, \textrm{s}^{-1}$ in the rest frame. No Doppler shift is applied to spectra during observations. This allows us to search for H i in galaxies with cz in the range of $-15\,908$ to $+126\,034\ \textrm{km} \, \textrm{s}^{-1}$ . Spectra can be recorded at maximum rate of 10 Hz, but we chose the lower rate of 1 Hz.

For each beam, three linear polarisation signals ( $\textrm{XX}^*$ , $\textrm{YY}^*$ and $\textrm{XY}^*$ ) are recorded in fits format files, each with a maximum size of 2 GB. No Doppler shift is applied to the raw spectra before the data is recorded. At this stage, the metadata includes the position of primary focus relative to the spherical centre, the skew of the lower platform, and the rotation angle of 19-beam receiver, which are all recorded in an independent spreadsheet. The sky coordinates and Doppler shifts are computed during the data reduction procedure by using the python package, ASTROPY (Astropy Collaboration et al. 2013, 2018).

3.2. Field selection

3.2.1. Target fields

In order to achieve our scientific goals, we apply the following criteria to select target fields:

1. 1. The Declination of the target fields should be close to $+25^{\circ}$ in order for FAST to utilise its full gain for the maximum length of observing time per day.

2. 2. The target fields should be evenly distributed in RA to minimise the influence of cosmic variance and maximise the ease of scheduling night observations at any time of year.

3. 3. Target fields should avoid strong radio continuum sources ( $S_{1.4}\textrm{GHz}>50$ mJy) to avoid reduction in sensitivity due to imperfect spectral subtraction. These are normally background sources, so there is no implicit selection bias.

4. 4. The target fields should overlap with SDSS and other imaging and spectroscopic surveys to provide corresponding optical identifications and galaxy parameters.

For the pilot FUDS0 field, we deliberately select a field which overlaps with the GAL2577 field of AUDS, which slightly compromises the continuum source criterion. In FUDS0, there are four continuum sources with flux densities larger than 50 mJy. For the remaining five target fields, all selection criteria are met. Currently, observations for two fields are in hand: the pilot FUDS0 field, and FUDS1. Details are listed in Table 1 and the field positions are illustrated in Figure 2.

Table 1. Positions of the FUDS target fields. Observations for FUDS0 are complete and observation for FUDS1 have commenced. Observations for FUDS2 – FUDS5 have not commenced. Their positions are preliminary.

Figure 2. This shows the six target fields in FUDS, shown by the black squares. The sky coverage of SDSS is shown with the grey area, while the DESI footprint (Levi et al. Reference Levi2019) is shown with the hatched area. The Galactic plane is shown in light blue.

3.2.2. Calibrators

Due to the deformable primary, the FAST gain barely changes at different zenith angles (within 26 $.^{\circ}$ 4) (Jiang et al. Reference Jiang2020). We only use one flux calibrator for each target field in FUDS using the following selection criteria:

1. 1. Strong ( $> 1.0\ \textrm{Jy}$ at 1.4 GHz) and stable flux density (marked as ‘P’ in the VLA calibrator list^{Footnote a} for the VLA D configuration, which has the beam size most similar to FAST).

2. 2. In order to avoid confusion, no other strong continuum sources close to the selected calibrator are listed in the NVSS catalogue (Condon et al. Reference Condon1998).

3. 3. Calibrators are close to target fields, and can be observed within a zenith angle of +26 $.\!\!^\circ$ 4.

Flux densities at frequencies near 1.0–1.5 GHz and their uncertainties were extracted from NED.^{Footnote b} To derive the frequency dependence of the calibrator flux densities, we use a power law to fit the measurements and uncertainties by using an optimised $\chi^2$ method:

(1) \begin{equation} S_{\nu} = S_{\nu_0} \left(\frac{\nu}{\nu_0}\right)^{\alpha}, \end{equation}

where the $\nu$ is frequency, $\nu_0$ is the reference frequency (1.4 GHz), and $\alpha$ is the spectral index. For FUDS, we assume the variability of compact sources is negligible during the course of our observations of a particular field.

As an example, we show the parameters for the FUDS0 calibrator, 0603+219, in Table 2 and show the fit to its continuum spectrum in Figure 3. The VLA flux density is 2.75 Jy at 1.4 GHz and the computed spectral index is $\alpha=0.565\pm 0.034$ . For computing frequency-dependent system parameters, we anchor the flux density with the accurate VLA measurement, and use our computed spectral index.

Figure 3. Flux densities and corresponding uncertainties for the calibrator at different frequency from previous works. The red line is the best fit of Equation (1). The 1- $\sigma$ uncertainty is shown by grey shadow.

Table 2. The properties of FUDS0 calibrator. The first four rows show the properties of the calibrator from VLA calibrator list. The flux density at different frequency is listed from the fifth to eleventh row. The parameters of Equation (1) from best fit are given in Row 12 and 13.

3.3. Observational strategy

3.3.1. Target fields

As shown in Figure 4 for the FUDS0 field, an on-the-fly (OTF) observation mode is employed. Each target field is $30'\times30'$ , as scanned by the central beam, with a rotation angle of $0^{\circ}$ . However, due to the multiple beam pattern, this results in a field size of $50'\times50'$ with good sensitivity. The spacing between scans is 1 arcmin, which is smaller than the Nyquist sampling interval for an individual 2.9 arcmin beam at 1.4 GHz. The target field is scanned in both RA and Dec direction with a rate of $15\ \textrm{arcsec} \, \textrm{s}^{-1}$ . Each scan takes 120 s and spectra are recorded every second. We use the high-temperature calibration noise diode (10 K) to inject noise for 1 s every 60 s. The total observing time per set of 31 scans is around 1.2 h (62 min scanning, 10 min overhead), or 2.4 h for the combined RA and Dec scans. Multiple sets of similar observations are then combined into the final data cube.

Figure 4. The position of the central beam whilst scanning the FUDS0 target field is shown by the red lines. The grey area has high sensitivity, but coverage extends to the outer black lines. The start positions for each of the 19 beams are given by the black circles.

4.1. RFI detection and baseline fitting

The GNSS signals are strong (having comparable intensity of the system temperature of FAST) and highly variable in time. They also occupy a large fraction of the FAST bandwidth (37.5%). The signals can affect the system temperature and its measurement even at non-RFI frequencies. The RFI from the compressor also adds uncertainty for flux intensity measurements by biasing the bandpass calibration. In order to mitigate these affect, we identify the RFI and to fit an underlying spectral baseline using the following steps:

1. 1. The spectrum is divided into several segments with an overlap of 15 MHz. The width of each segment varies from 90 to 240 MHz, depending on the number of channels that are free of RFI.

2. 2. For each segment, a Gaussian weighted polynomial function is employed to fit the baseline. The standard deviation of the Gaussian weight is 0.17 times the segment width. The order of the fit is iterated from 0 to 6. For each order, a pseudo-RMS is calculated using 1.4826 $\times$ MAD (Median Absolute Deviation). The channels with absolute values exceeding $7\times$ RMS and immediate neighbouring channels with absolute values exceeding $3\times$ RMS are marked as RFI.

3. 3. A temporary baseline for the whole spectrum is determined by combining the baselines of each segment with Gaussian weights.

4. 4. The above procedures are iterated on a smoothed baseline-removed spectrum with increasing smoothing length (25, 50, 100 and 200 kHz). The intermediate baseline is generated by summing together the temporary baselines from each iteration.

5. 5. The final baseline is the intermediate baseline generated in the final iteration. All the channels ever marked as RFI in each iteration are recorded as RFI in the final spectrum.

We note that RFI flagging in raw data is time consuming. Since the RFI is stable in frequency in a short period (around 20 min, see Section 4.2), we only perform this procedure on one spectrum per 60 spectra (the same frequency as the noise injection in the target fields). The template for flagged RFI is therefore used for the next 59 spectra.

In Figure 5, we show a total power spectrum ( $\textrm{XX}^*$ polarisation) from a single 1 s integration from the FUDS0 central beam. The fitted baseline and the identified regions with RFI are demonstrated. About 41.7% of the spectrum contains RFI at some level. This procedure may miss weak RFI, but the baseline fitting is relatively robust.

Figure 5. The upper panel shows a total power spectrum (polarisation $\textrm{XX}^*$ ) from the central beam (black line) and the underlying fitted baseline (red line). The units are arbitrary. The spectrum around 1.23 GHz is zoomed in to show the standing wave. The frequencies of the detected RFI are indicated by the blue bars in the lower panel.

4.2. Compressor RFI

To characterise the signals generated by the refrigerating dewar in the compressor, we analyse the calibrator data taken on 2019 August 25. In Figure 6, we show a waterfall image of the calibrated spectra between frequencies of 1.31 and 1.41 GHz. This frequency range excludes the GNSS RFI and H i emission from the Milky Way. The compressor signals (the vertical stripes) generally have an intensity of about 0.3 Jy and a linewidth of ${\sim}1$ MHz. The signals are variable in amplitude and frequency (middle panel of Figure 6) and also shift in frequency when comparing data on different days (see Figure 7). The signals do not depend on azimuth or elevation as the compressor is housed in the focus cabin.

Figure 6. In the left panel, calibrated data between 1.31 and 1.41 GHz is shown for the $\textrm{XX}^*$ polarisation channel of beam 1 as a function of frequency (in GHz) and time (in sec). The right panel is a zoom (in frequency space) of the region indicated by the red rectangle. The signals present are due to RFI from a compressor in FAST data taken prior to 2021 July.

Figure 7. The spectra of one polarisation ( $\textrm{XX}^*$ ) from beam 1 on 2019 Aug 25 (red) and 2019 Dec 28 (blue). The former spectrum has been vertically shifted by $0.5\ \textrm{Jy beam}^{-1}$ for clarity.

Calibrated spectra of all 19 beams from the same observation are shown in the left panel of Figure 8. There is a similarity in frequency and intensity of the compressor RFI across all 19 beams. We create a smoothed median spectrum from all beams and use it as template to subtract from all spectra. The residual spectra after RFI removal are shown in the right panel of Figure 8. Removal of the compressor signals is only performed after the calibration due to large power differences between the beams. Large sources which extend across several beams will potentially have their flux density reduced. We note that further flagging is required to remove residual RFI.

Figure 8. Example spectra of all beams observed at one time (polarisation $\textrm{XX}^*$ from beams 1 to 19 are offset from the bottom to the top) are shown in the left panel. The frequency and intensity of the local RFI is fairly constant across all beams. The right panel shows the residual spectra after removing a template constructed from all beams. The RFI intensity has been reduced.

5.1. Gain stability

The FAST calibration method involves switching on and off a noise diode, which is assumed to have a constant amplitude. Once the spectrum of noise diode has been calibrated against a radio continuum source, the diode noise can then be used to measure the system equivalent flux density (SEFD) in Jy, $S_\textrm{sys}$ and its stability. Because the high-temperature noise diode elevates $T_{\rm sys}$ by nearly 50%, it is impractical to use with a 50% duty cycle for target field observations. We therefore needed to confirm that the receiver calibration was stable in a mode where the calibration ‘on’ spectrum was only taken for 1 s every 60 s.

To confirm this, we first removed the identified RFI channels (see Section 4.1), and channels that are outside of the official bandpass ( $< 1.05$ , or $> 1.45$ GHz). The noise diode was switched on and off at a rate of 1 Hz while performing a calibrator observation, and the difference in the mean power in each pair of on/off spectra is calculated. Figure 9 shows the result, which indicates that no significant variation occurs over any time frame less than 1 500 s. This confirms that both the noise diode and the receiver gain are adequately constant within this period.

Figure 9. The power difference between successive integration cycles with the calibration noise diode on and off, respectively, as a function of time. The power is calculated after removing edge channels and channels affected by RFI. The units of power are arbitrary. The grey area indicates a time range where the calibrator distance was less than 3 arcmin.

Figure 10. The mean power (in arbitrary units) of polarisation $\textrm{XX}^*$ from beam 1 as a function of angular distance from the calibrator. Channels containing RFI and edge channels are removed prior to calculating the mean. The mean power for integration cycles with the noise diode on and off are shown with red plus signs and blue dots, respectively. The best fits curves using Equation (2) are shown.

5.2. Instrumental characterisation

For a radio telescope, a Gaussian function is a good approximation for angular response of the beam (the point spread function). The following function can therefore be used to fit a 1-dimensional scan of a calibrator where the noise diode (CAL) is either on or off:

(2) \begin{equation} f(d)= \left\{ \begin{array}{ll} Ae^{-\dfrac{(d-d_o)^2}{2\sigma^2}}+C_1 &\text{(CAL on)}\\ \\[-7pt] Ae^{-\dfrac{(d-d_o)^2}{2\sigma^2}}+C_2 &\text{(CAL off)} \end{array} \right. \end{equation}

where d is the angular distance offset from the calibrator, $d_o$ is the calibrator position, A is the amplitude, $\sigma$ is the standard deviation, and $C_1$ and $C_2$ are the CAL on and off baselines, respectively. The flux density of the noise diode is then $S_\textrm{noise}=\frac{C_1-C_2}{A}\times S_\textrm{calibrator}$ , the SEFD $S_\textrm{sys}=\frac{C_2}{A}\times S_\textrm{calibrator}$ , and the full width half power (FWHP) beamwidth is $\theta=\sqrt{8\ln2} \sigma$ , where $S_\textrm{calibrator}$ is the flux density of the calibrator. We also measure the pointing error, $\sigma_\textrm{pointing}$ from the RMS of the deviation of $d_o$ from the catalogued position.

We firstly calculate mean instrumental parameters, $\bar{S}_\textrm{noise}$ , $\bar{S}_\textrm{sys}$ , $\bar{\sigma}_\textrm{pointing}$ , and $\bar{\theta}$ averaged over the frequency range 1.05–1.45 GHz, using the VLA calibrator flux density at 1.4 GHz.

However, we also derive frequency-dependent parameters, $S_\textrm{noise}(\nu)$ , $S_\textrm{sys}(\nu)$ , $\sigma_\textrm{pointing}(\nu)$ , and $\theta(\nu)$ using the frequency-dependent flux density of calibrator, $S_\textrm{calibrator}(\nu)$ derived in Section 3.2.2. Some uncertainty is introduced at frequencies affected by RFI, so stable results can only be derived using the baseline of raw spectrum provided in Section 4.1. In some cases (e.g. for beam size), the frequency dependence is strong, as expected. In other cases (e.g. pointing error), there is virtually no frequency dependence. An example Gaussian fit to the beam response for polarisation $\textrm{XX}^*$ , beam 1 is shown in Figure 10, and the derived instrumental parameters are shown in Figure 11. Due to strong GNSS RFI the measurements between 1.15 and 1.30 GHz (the grey area Figure 11) are less accurate.

Figure 11. The variation of noise diode flux density (top left), SEFD (top right), pointing error (bottom left) and beam size (bottom right) as a function of frequency is shown in black for beam 1, polarisation $\textrm{XX}^*$ . The hatched areas show the 1- $\sigma$ uncertainties. The grey area shows the frequency range from 1.15 to 1.3 GHz which is affected by RFI (GNSS satellites). The mean values, which exclude RFI-affected regions and edge channels, are shown in red.

5.2.2. System equivalent flux density (SEFD)

In the upper right panel of Figure 11, we show the SEFD. The mean value is $\bar{S}_\textrm{sys}=1.24 \pm 0.08\ \textrm{Jy}$ , but $S_\textrm{sys}(\nu)$ increases by 15% below 1.05 GHz, which is outside of the nominal frequency range. A bump is also seen between 1.08 and 1.12 GHz. Similar deviations are found in other beams and in the $\textrm{YY}^*$ polarisation. Since it also seems to correspond to a change in the SEFD and beamwidth, it is possibly an instrumental effect in the feed or the feed-LNA transition.

In Figure 12, we derive a median frequency-dependent correction factor for the mean SEFD correction by using data from both polarisation on all 19 beams. The correction factor is parameterised according to:

(3) \begin{equation} f_{\rm sys}(\nu) = A \times e^{-(\frac{\nu - \nu_0}{B})} + C \times e^{-(\frac{\nu - \nu_1}{D})^2} + E, \end{equation}

where the first term is the exponential rise at low frequency ( $< 1.05$ GHz), the second term is the Gaussian component at 1.1 GHz, and the third term is the deviation at high frequency ( $> 1.3$ GHz). We use the optimised $\chi^2$ method to derive the best fit parameters for the median line, excluding the RFI affected frequency range (1.15–1.30 GHz). The derived best fit parameters are given in Table 3 and the best fit line is shown by the red solid line. The uncertainty for each of the parameters is estimated using a jackknife method, that is removing one spectrum from 19 beams $\times$ 2 polarisation channels each time.

Figure 12. The median normalised SEFD across all beams and polarisation channels is shown as a function of frequency (black dashed line), with corresponding 1- $\sigma$ uncertainty (hatched area). A rise is apparent at low frequency ( $< 1.15$ GHz) in all beams. We use Equation (3) to fit the median, and derive the best fit line shown by the solid red line. Data affected by RFI (grey area) are excluded in the fit.

Table 3. The best fit parameters for the normalised frequency dependence of the SEFD using Equation (3). The uncertainty is from the jackknife method.

To evaluate the effects of solar interference on the SEFD, we compare the $\bar{S}_{\textrm{sys}}$ of each beam for the day and night observations (see Figure 13). We find values for $\bar{S}_{\textrm{sys}}$ which are about 16% higher in the day. The mean $\bar{S}_{\textrm{sys}}$ values are 1.34 and 1.15 Jy in the day and night, respectively.

Figure 13. The SEFD for all 19 beams is shown. The red pluses represent day time data, while the blue dots represent night time data. The beams in outer circles of the receiver have higher SEFD. The mean daytime SEFD is 1.34 Jy, and 1.15 Jy at night.

For daytime observations, we list the mean $\bar{S}_{\textrm{sys}}$ and its uncertainty for different combinations of beams in Table 4. The central beam (beam 1) has the lowest $\bar{S}_{\textrm{sys}}$ value of 1.23 Jy. The inner circle of beams (beams 2–7) has a similar mean value ( $S_{\rm sys} = 1.29\ \textrm{Jy}$ ). The second circle (beams 9, 11, 13, 15, 17 and 19) has a significantly higher mean value ( $\bar{S}_{\rm sys}=$ 1.36 Jy), and the final circle (beams 8, 10, 12, 14, 16 and 18) have the highest mean value ( $\bar{S}_{\rm sys} = 1.40\ \textrm{Jy}$ ). Table 4 shows that this is not due to a rise in system temperature $T_{\textrm{sys}}$ , which is fairly constant (here we adopt the gain at 1.35 GHz cited by Jiang et al. Reference Jiang2020 to convert SEFD to $T_{\textrm{sys}}$ ). Instead, it is due to the classical drop in off-axis antenna efficiency for paraboloid antennas (Ruze Reference Ruze1965).

Table 4. The mean SEFD, system temperature and mean beam size for beams at different offsets from the optical axis (beam 1) for observations of 2019 Aug 25.

5.2.3. Pointing error

In the lower left panel of Figure 11, we show the measured pointing error. It clearly shows that $\sigma_\textrm{pointing}(\nu)$ is nearly constant, and consistent with the mean value, $\bar{\sigma}_\textrm{pointing}$ . We also show the histogram of mean pointing errors across all beams and polarisations over two days of observations (76 measurements) in Figure 14. The mean and median errors are 6.3 and 8.1 arcsec, respectively. Considering the 2.8 arcmin beam size of FAST at 1.5 GHz, the pointing accuracy is better than $\frac{1}{20}$ of the FAST beam.

Figure 14. A histogram of the pointing errors measured in both polarisations, all beams, and in both observing days. The bins have a width of 5 arcsec. The median and mean values are 6.30 and 8.05 arcsec, and are shown by red dotted line and red dashed line, respectively.

5.2.4. Beam size

The lower right panel of Figure 11 shows the beam size measurements. The mean beam size $\bar{\theta}$ , is 3.26 arcmin. The beam size for beams at different radii from the optical axis show a weakly increasing trend (see the final column in Table 4). Our findings are consistent with the results from Jiang et al. (Reference Jiang2020).

As expected, the beam size $\theta(\nu)$ , decreases with frequency. Across all beams and both polarisations, the median frequency-dependent beam size and its 1- $\sigma$ uncertainty is shown in Figure 15. The beam size is dependent on wavelength and telescope diameter. For an unblocked circular aperture with uniform illumination, the FWHP beam size is $\theta = 1.02 \lambda/D$ in radians (Wilson, Rohlfs, & Hüttemeister Reference Wilson, Rohlfs and Hüttemeister2013), where $\lambda$ is the wavelength and D is the illuminated diameter. Normal radio telescopes with tapered illumination patterns will typically have beam sizes which are larger by 20%, but this depends on the taper (Wilson et al. Reference Wilson, Rohlfs and Hüttemeister2013) and the offset from the optical axis. For the similar prime-focus Parkes multibeam system (Staveley-Smith et al. Reference Staveley-Smith1996), the average FWHP beam size is $\theta = 1.225 \lambda/D$ in radians. Hence, for the multiple beams on FAST, we parameterise the beam size as follows:

(4) \begin{equation} \theta = 1.225 \frac{\lambda}{A + B\left( \frac{\lambda-0.21}{0.21}\right) + C\left(\frac{\lambda-0.21}{0.21}\right)^2}, \end{equation}

where A, B and C are functions of $\lambda$ . We use an optimised $\chi^2$ method to find the best fit (see Figure 15). The best fit parameters and the parameters for all the 19 beams are listed in Tables 5 and Tables A.1, respectively.

Figure 15. The median FWHP beam size for all 19 beams is shown by the black solid line, with the corresponding 1- $\sigma$ uncertainty shown by the hatched area. The best fit line by using Equation (4) is shown by the red dashed line. The beam size in our final cube (post-gridding, described in Section 6.4) is shown by the red solid line.

Table 5. Parameters of Equation (4) for the beam size of FAST, and beam size in our final cube.

6.1. System equivalent flux density

For target field flux density calibration, we use the flux density of noise diode as calibrated by the above calibrator observations. We first remove channels affected by the known RFI and channels outside the FAST nominal bandpass ( $< 1.05$ or $> 1.45$ GHz) (for calibration purposes only). Subsequently, we calculate the mean power for all the spectra. We derive the CAL step by subtracting the CAL off spectrum from the CAL on spectrum that it follows. The value for $\bar{S}_\textrm{noise}$ derived in Section 5.2 is used as the equivalent flux density of the CAL step. The derived value for $\bar{S}_{\textrm{sys}}$ is then used for the subsequent 59 spectra. The frequency-dependent scale factor derived in Section 5.2 is applied after flagging.

6.2. Bandpass removal

We follow the robust bandpass calibration procedure as used for HIPASS (Barnes et al. Reference Barnes2001). In order to calibrate a raw spectrum $P_i(\nu)$ according to its measured SEFD $S_i$ , we select neighbouring spectra as a reference spectrum set, $C_r$ . The reference SEFD $S_{\textrm{sys},r}$ , and the reference spectrum $P_r(\nu)$ , are computed from the median of the spectra in $C_r$ . Subsequently, the raw spectrum is flattened using the reference spectrum and scaled to units of $S_{\textrm{sys},r}$ . The system temperature contribution, including sky continuum, is then removed by subtracting $S_{sys, i}$ . The procedure is summarised by:

(5) \begin{equation} S_i(\nu) = S_{\textrm{sys}, r} \times \frac{P_i(\nu)}{P_r(\nu)} - S_{\textrm{sys}, i}, \end{equation}

where $S_i(\nu)$ is the calibrated spectrum in Jy, $S_{\textrm{sys}, i}$ and $S_{\textrm{sys}, r}$ are in Jy, and $P_i(\nu)$ and $P_r(\nu)$ are in (arbitrary) units of power. For each beam, we include the nearest 260 spectra in the same beam as the reference spectrum set, giving spatial coverage of about two scan lines. As explained in Calabretta , Staveley-Smith, & Barnes (Reference Calabretta, Staveley-Smith and Barnes2014), the form of Equation (5) flattens the bandpass, converts the flux scale to Jy and removes continuum to zeroth-order. Due to residual spectral artefacts (spectral index and ripple) in strong continuum sources ( $> 10$ mJy), spectra within 3 arcmin of these sources are excluded in the reference spectrum set. To reduce sidelobes caused by any strong HI emission, we use a 5-channel Hanning window (0.25, 0.75, 1, 0.75, 0.25) to convolve with the calibrated spectra. Therefore, the frequency resolution is lowered to 22.9 kHz. This procedure is performed for each polarisation channel and for all 19 beams. We use a cubic polynomial function to fit the calibrated spectra to further flatten the baseline. Figure 16 shows 60 calibrated spectra for the $\textrm{XX}^*$ polarisation of beam 1. The residual compressor RFI, particularly the vertical strips between 1.31 and 1.41 GHz, are still visible (upper left panel). The use of the median spectrum for all 19 beams suppresses the compressor RFI (upper right panel) but also the extended Galactic H i emission.^{Footnote c}

Figure 16. Upper left panel: calibrated data for polarisation $\textrm{XX}^*$ , beam 1 in the frequency-time domain. Upper right panel: constant RFI is suppressed by removing the median spectrum across the 60 integration cycles. Lower left: data flagged as RFI is removed. Lower right panel: final frequency-dependent intensity calibration is applied.

6.3. RFI removal

To remove residual RFI from all target spectra, we adopt the following flagging steps:

1. 1. Pixels without extreme values (not deviating more than $\pm$ 0.1 Jy) are used for deriving noise. The theoretical sensitivity for the given frequency and time resolution is 15.3 mJy, assuming a SEFD of 1.3 Jy (corresponding to a gain of $15\ \textrm{K}\ \textrm{Jy}^{-1}$ and a system temperature of 20 K).

2. 2. Channels with high noise levels, $\sigma_c$ > $1.2\sigma_g$ are noted as bad, where $\sigma_c$ is computed as 1.4826 $\times$ MAD for each channel in all scans in one dimension (either RA or Dec) over one day, and $\sigma_g$ is the global average calculated in the same way over all channels.

3. 3. Pixels in bad channels with high noise levels in the time-frequency domain are flagged, that is $\sigma_p$ > $2\sigma_g$ , where $\sigma_p$ is defined as the RMS in the surrounding $5 \times 5$ box. This step also detects regions where the mean flux density deviates from zero, so bright galaxies in RFI-affected channels may be eliminated.

4. 4. Any flagged pixel is grown into a flagged $5 \times 5$ pixel region.

5. 5. The spectra take with CAL on, which have higher $S_{\textrm{sys}}$ , are also flagged.

The third panel of Figure 16 shows the effect of removing the flagged pixels. The flagged data accounts for 25.6%. of the calibrated data, which is smaller than the 41.7% of flagged RFI data in the raw spectral data (Figure 5). This noise-based flagging method therefore recovers 15% of the data that would otherwise have been flagged. This is partly helped by the compressor RFI template. As already noted, the frequency-dependent flux density correction factor, $f_\textrm{sys}(\nu)$ is applied to the data after flagging, and also shown in Figure 16.

6.4. Gridding

We employ the robust gridding algorithm developed for HIPASS to create our final FUDS data cube. The final cube has a size of $1^{\circ} \times 1^{\circ}$ and a pixel size of 1 arcmin $\times$ 1 arcmin. The frequency spacing in the final data cube is the same as the raw spectra, 7.63 kHz, but the resolution after Hanning smoothing is 22.9 kHz. For each pixel, we use a radius of 1.3 arcmin as a cutoff to select the contributing spectra. Considering the variation of beam size with frequency, up to 40% across the whole bandwidth, the median frequency-dependent beam size derived in Section 5.2 is used for all beams. The spectrum for a given position in the cube is then derived as follows:

(6) \begin{equation} S(\nu)=\frac{\textrm{median}(S_r(\nu))}{\textrm{median}(w_r(\nu))}, \end{equation}

where $S_r(\nu)$ are the spectra within 1.3 arcmin and $w_r(\nu)$ are their corresponding Gaussian weights.

Assuming all the 19 beams have same analytic beam size, $\theta(\nu)$ , derived in Section 5.2, we perform the same gridding procedure above to derive the beam size in our final cube by using a Monte Carlo (MC) method. We then use Equation (4) to fit the gridded beam. The best fit line is shown in Figure 15, and the best fit parameters of Equation (4) are given in Table 5. This beam size is required when deriving precise flux densities for extended source in the gridded cubes.

6.5. Strong continuum sources

We employ the method of Barnes et al. (Reference Barnes2001) to remove residual emission near strong continuum sources. We use the spectra toward continuum sources with flux densities stronger than 10 mJy to generate a template spectrum, given by:

(7) \begin{equation} \textbf{S}_T(\nu) = \Sigma_i w_i S_i(\nu), \end{equation}

where $S_i(\nu)$ is the central spectrum of ith continuum source, and $w_i$ is the weight. The spectrum associated with the strongest continuum source has a weight set to unity. The weight of other sources is determined by the slope of the linear regression between two spectra $S_i(\nu)$ and $S_1(\nu)$ in the frequency ranges 1.015–1.137 GHz, 1.30–1.416 GHz and 1.423–1.463 GHz. For each spectrum in the final data cube, we perform a linear regression with the template in the same manner to derive the weight w, then subtract $w\textbf{S}_T(\nu)$ from the spectrum.