<figure><img src="image_1.png" /></figure>

HAL
open science

# FAST drift scan survey for HI intensity mapping: I. preliminary data analysis

Yichao Li, Yougang Wang, Furen Deng, Wenxiu Yang, Wenkai Hu, Diyang Liu, Xinyang Zhao, Shifan Zuo, Shuanghao Shu, Jixia Li, et al.

► To cite this version:

Yichao Li, Yougang Wang, Furen Deng, Wenxiu Yang, Wenkai Hu, et al. FAST drift scan survey for HI intensity mapping: I. preliminary data analysis. The Astrophysical Journal, 2023, 954 (2), pp.139. 10.3847/1538-4357/ace896. hal-04103969

HAL Id: hal-04103969

https://hal.science/hal-04103969v1

Submitted on 27 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# FAST Drift Scan Survey for HI Intensity Mapping: I. Preliminary Data Analysis

Yichao Li<fnref n="1" />, Yougang Wang<fnref n="2" />,<fnref n="3" />,<fnref n="4" />, Furen Deng<fnref n="2" />,<fnref n="3" />,<fnref n="4" />, Wenxiu Yang<fnref n="2" />,<fnref n="3" />,<fnref n="4" />, Wenkai Hu<fnref n="5" />,<fnref n="6" />, Diyang Liu<fnref n="1" />, Xinyang Zhao<fnref n="1" />, Shifan Zuo<fnref n="2" />,<fnref n="4" />, Shuanghao Shu<fnref n="2" />,<fnref n="3" />, Jixia Li<fnref n="2" />,<fnref n="4" />, Peter Timbie<fnref n="7" />, Réza Ansari<fnref n="8" />, Olivier Perdereau<fnref n="9" />, Albert Stebbins<fnref n="10" />, Laura Wolz<fnref n="11" />, Fengquan Wu<fnref n="2" />,<fnref n="3" />,<fnref n="4" />, Xin Zhang<fnref n="1" />,<fnref n="12" />,<fnref n="13" />, and Xuelei Chen<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" />

Received 2023 May 10; revised 2023 July 7; accepted 2023 July 17; published 2023 August 31

## Abstract

This work presents the initial results of the drift-scan observation for the neutral hydrogen (HI) intensity mapping survey with the Five-hundred-meter Aperture Spherical radio Telescope (FAST). The data analyzed in this work were collected in night observations from 2019 through 2021. The primary findings are based on 28 hr of drift-scan observation carried out over 7 nights in 2021, which covers 60 deg<fnref n="2" /> sky area. Our main findings are, first, our calibration strategy can successfully correct both the temporal and bandpass gain variation over the 4 hr drift-scan observation. Second, the continuum maps of the surveyed region are made with frequency resolution of 28 kHz and pixel area of 2.95 arcmin<fnref n="2" />. The pixel noise levels of the continuum maps are slightly higher than the forecast assuming Tsys = 20 K, which are 36.0 mK (for 10.0 s integration time) at the 1050–1150 MHz band, and 25.9 mK (for 16.7 s integration time) at the 1323–1450 MHz band, respectively. Third, the flux-weighted differential number count is consistent with the NRAO-VLA Sky Survey (NVSS) catalog down to the confusion limit ~7 mJy beam⁻<fnref n="1" />. Finally, the continuum flux measurements of the sources are consistent with those found in the literature. The difference in the flux measurement of 81 isolated NVSS sources is about 6.3%. Our research offers a systematic analysis for the FAST HI intensity mapping drift-scan survey and serves as a helpful resource for further cosmology and associated galaxies sciences with the FAST drift-scan survey.

*Unified Astronomy Thesaurus concepts: Large-scale structure of the universe (902); Astronomy data analysis (1858); Surveys (1671)*

## 1. Introduction

Measurements of the cosmological large-scale structure (LSS) play an important role in studying the evolution of the Universe. In the past decades, the LSS fluctuations have been explored by observing the galaxy distribution in the Universe with wide-field spectroscopic and photometric surveys (Cole et al. 2005; Eisenstein et al. 2005; Anderson et al. 2014; Hinton et al. 2017; Alam et al. 2021). Recently, another cosmological probe of the LSS by observing the neutral hydrogen (HI) in the galaxies via its 21 cm emission line of hyperfine spin-flip transition has been proposed (e.g., Battye et al. 2004; McQuinn et al. 2006; Pritchard & Loeb 2012).

Several HI galaxy surveys have already been carried out, e.g., the 64 m Parkes telescope in Australia with the HI Parkes All-Sky Survey (HIPASS; Barnes et al. 2001; Meyer et al. 2004;

Zwaan et al. 2004), the 76 m Lovell Telescope at Jodrell Bank with the HI Jodrell All-Sky Survey (HIJASS; Lang et al. 2003), the Arecibo Legacy Fast ALFA (ALFALFA) survey (Giovanelli et al. 2005, 2007; Saintonge 2007), and the Jansky Very Large Array (JVLA) 10 deg<fnref n="2" /> deep survey (Jarvis et al. 2014). However, limited by the sensitivity and the angular resolution of the radio telescopes, the redshift range of these surveys is much smaller than the current optical surveys. To resolve the HI emission line from individual distant galaxies at centimeter wavelength requires a large radio interferometer and is also time-consuming. Instead, a technique known as HI intensity mapping (HI IM), which measures the total HI intensity of many galaxies within large voxels (Chang et al. 2008; Loeb & Wyithe 2008; Mao et al. 2008; Pritchard & Loeb 2008; Wyithe & Loeb 2008; Wyithe et al. 2008; Peterson et al. 2009; Bagla et al. 2010; Seo et al. 2010; Lidz et al. 2011; Ansari et al. 2012; Battye et al. 2013), can be quickly carried out and extended to very large survey volume, and is ideal for cosmological surveys (Xu et al. 2015; Jin et al. 2021; Zhang et al. 2021; Wu et al. 2023a, 2023b; Wu & Zhang 2022; Zhang et al. 2023).

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Key Laboratory of Cosmology and Astrophysics (Liaoning) & College of Sciences, Northeastern University, Shenyang 110819, People's Republic of China; wangyg@bao.ac.cn, zhangxin@mail.neu.edu.cn</footnote>

<footnote n="2">² National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, People's Republic of China; xuelei@cosmology.bao.ac.cn</footnote>

<footnote n="3">³ School of Astronomy and Space Science, University of Chinese Academy of Sciences, Beijing 100049, People's Republic of China</footnote>

<footnote n="4">⁴ Key Laboratory of Radio Astronomy and Technology, Chinese Academy of Sciences, A20 Datun Road, Chaoyang District, Beijing 100101, People's Republic of China</footnote>

<footnote n="5">⁵ Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="6">⁶ ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia</footnote>

<footnote n="7">⁷ Department of Physics, University of Wisconsin—Madison, Madison, Wisconsin 53706, USA</footnote>

<footnote n="8">⁸ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, F-91191, Gif-sur-Yvette, France</footnote>

<footnote n="9">⁹ Université Paris-Saclay, CNRS/IN2P3, IJCLab, F-91405 Orsay, France</footnote>

<footnote n="10">¹⁰ Fermi National Accelerator Laboratory, P.O. Box 500, Batavia IL 60510, USA</footnote>

<footnote n="11">¹¹ Jodrell Bank Centre for Astrophysics, Department of Physics & Astronomy, The University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="12">¹² National Frontiers Science Center for Industrial Intelligence and Systems Optimization, Northeastern University, Shenyang 110819, People's Republic of China</footnote>

<footnote n="13">¹³ Key Laboratory of Data Analytics and Optimization for Smart Industry (Ministry of Education), Northeastern University, Shenyang 110819, People's Republic of China</footnote>
The HI IM technique has been explored by measuring the cross-correlation function between an HI IM survey carried out with Green Bank Telescope (GBT) and an optical galaxy survey (Chang et al. 2010). Later, a few detections of the cross-correlation power spectrum between an HI IM survey and an optical galaxy survey were reported with the GBT and the Parkes telescope (Masui et al. 2013; Wolz et al. 2017; Anderson et al. 2018; Wolz et al. 2022; Amiri et al. 2023). Several ongoing HI IM experiments focus on the post-reionization epoch, such as the Tianlai project (Chen 2012; Li et al. 2020; Wu et al. 2021; Perdereau et al. 2022; Sun et al. 2022) and the Canadian Hydrogen Intensity Mapping Experiment (CHIME, Bandura et al. 2014). Several HI IM experiments are under construction, such as the Baryonic Acoustic Oscillations from Integrated Neutral Gas Observations (BINGO, Battye et al. 2013) and the Hydrogen Intensity and Real-Time Analysis experiment (HIRAX, Newburgh et al. 2016). The HI IM technique is also proposed as the major cosmology project with the Square Kilometre Array (SKA;<fnref n="14" /> Santos et al. 2015; Square Kilometre Array Cosmology Science Working Group et al. 2020) and MeerKAT (Bull et al. 2015; Santos et al. 2016; Li et al. 2021; Paul et al. 2021; Wang et al. 2021; Chen et al. 2023). Recently, the MeerKAT HI IM survey reported cross-correlation power spectrum detection with the optical galaxy survey (Cunnington et al. 2023). Meanwhile, using the MeerKAT interferometric observations, Paul et al. (2023) reports the HI IM auto power spectrum detection on Mpc scales. However, the HI IM auto power spectrum on large scales remains undetected (Switzer et al. 2013).

The Five-hundred-meter Aperture Spherical radio Telescope (FAST, Nan et al. 2011; Li & Pan 2016) was recently built and is now operating for observations (Jiang et al. 2020). FAST is located in the Dawodang karst depression, a natural basin in Guizhou province, China (E 106°86, N 25°65). FAST can reach all parts of the sky within 40° from the zenith, corresponding to ~25000 deg² sky area. With a 300 m effective aperture diameter, FAST is the most sensitive radio telescope in the world.

The observation efficiency of FAST is highly increased using the feed array, which is mounted in the focus cabin that is suspended on cables above the reflector. The feed array has 19 hexagonally arranged and closely spaced feeds, which point to different directions in the sky. The relative positions of the 19 feeds are shown in Figure 1. Currently, the receiver is working at the L-band, i.e., frequencies between 1.05 and 1.45 GHz. The feed array increases the field of view, which is ideal for a large-area survey. Additionally, with ~3' angular resolution, the large-area survey with FAST can potentially resolve the HI environment in a large number of galaxy clusters, groups, filaments, and voids.

Forecasts and simulations show that the large-area HI IM survey with FAST is ideal for cosmological studies (Li & Ma 2017; Hu et al. 2020). Before committing to a large survey, we proposed a pilot survey. With the pilot survey, we aim to find the relevant characteristics of FAST receiver system. The systematic 1/f-type gain variation is studied in Hu et al. (2021). In this work, we address the time-ordered data (TOD) analysis pipeline for the HI IM with FAST drift-scan observation. This paper is organized as follows. We describe the data collected

<figure><img src="image_5.png" /><figcaption>**Figure 1.** The positions of the feeds in the FAST L-band 19-feed array. The feed array without rotation is shown in the left-hand panel and the feed array rotated by 23°4 is shown in the right-hand panel. The circles represent the 3' beam size at a frequency of 1420 MHz. The central, inner-circle, and outer-circle feeds are shown in red, blue, and green, respectively.</figcaption></figure>

for HI IM pilot survey in Section 2 and the TOD analysis method in Section 3. The results and implications of the TOD analysis are discussed in Section 4. In Section 5, a summary of this work is presented.

## 2. Observational Data

The data analyzed in this work were collected in night sessions spanning 2019, 2020, and 2021. During the observation of each night, the telescope was pointed at a fixed altitude angle. We chose an area close to the zenith at the telescope site, corresponding to decl. ~26°. The resulting decl. in the J2000 equatorial coordinate varies slightly (at a magnitude of ~3 arcmin) during the 4 hr drift-scan.

During the observations in 2019 and 2020, we investigated different configurations of the observation parameters, e.g., the frequency resolution, sampling rate, feed array rotation angle and the level of the noise diode power, to test the impacts on the observation results. The finer resolution in frequency and time is beneficial for identifying narrow radio frequency interference (RFI) and resolving the spectroscopic profile of HI in galaxies. However, it produces a significant amount of data and takes a long time to process. Our testing results show that, with the exception of a few badly contaminated frequency ranges, RFI contamination can be properly detected utilizing the frequency resolution of 7.6 kHz. Meanwhile, the spectral features of interest can be resolved using this frequency resolution. Data with finer frequency resolution and higher sampling rates are used to measure the systematic 1/f noise at different temporal and spectroscopic scales. Our analysis (Hu et al. 2021) shows that the systematic 1/f noise is negligible within a few hundred seconds timescale. Thus, we adopt 7.6 kHz frequency resolution and 1 s integration time for the survey observations.

We use the FAST L-band 19-feed receiver. The position of the feed in the 19-feed array is shown in Figure 1. The feed array is rotated by 23°4 to obtain the maximum span of decl. coverage during the drift scans (Li et al. 2018). We also have a couple of observations without rotating the feed array. During such drift-scan observations, the same sky stripe is repeatedly scanned by different feeds, which is ideal for systematic checking via cross-correlating the observation of different feeds.

A noise diode is built into the receiver system and its output can be injected as a real-time calibrator during the observation. In every 8 s, the noise diode was fired for 0.9 s, which is

<footnote n="14">¹⁴ https://www.skao.int</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** The footprint of the FAST HI IM Pilot Survey observation using Zenith equal area projection (ZEA). The filled blue region shows the sky area of the major data observed in 2021, the green stripe indicates the test observation carried out on 2019 May 31st and 2020 May 8th, and the orange region shows the full target area for FAST HI IM Pilot Survey.</figcaption></figure>

<table><caption>**Table 1**
Observation Mode</caption><thead><tr><th>Field Center</th><th>Date</th><th>Frequency Resolution [kHz]</th><th>Integration Time [s]</th><th>Noise Diode Level</th><th>Rotation Angle [°]</th></tr></thead><tbody><tr><td>HIMGS 1100+2539</td><td>2019-05-27</td><td>0.5</td><td>1</td><td>high</td><td>0</td></tr><tr><td>HIMGS 1100+2554</td><td>2019-05-28</td><td>0.5</td><td>0.1</td><td>high</td><td>0</td></tr><tr><td>HIMGS 1100+2609</td><td>2019-05-29</td><td>7.6</td><td>0.1</td><td>low</td><td>0</td></tr><tr><td>HIMGS 1100+2554</td><td>2019-05-30</td><td>7.6</td><td>1</td><td>low</td><td>0</td></tr><tr><td>HIMGS 1100+2639</td><td>2019-05-31</td><td>7.6</td><td>1</td><td>low</td><td>23.4</td></tr><tr><td>HIMGS 1100+2639</td><td>2020-05-08</td><td>7.6</td><td>1</td><td>low</td><td>23.4</td></tr><tr><td>HIMGS 1100+2600</td><td>2021-03-02</td><td>7.6</td><td>1</td><td>low</td><td>23.4</td></tr><tr><td>HIMGS 1100+2632</td><td>2021-03-05</td><td>7.6</td><td>1</td><td>low</td><td>23.4</td></tr><tr><td>HIMGS 1100+2643</td><td>2021-03-06</td><td>7.6</td><td>1</td><td>low</td><td>23.4</td></tr><tr><td>HIMGS 1100+2654</td><td>2021-03-07</td><td>7.6</td><td>1</td><td>low</td><td>23.4</td></tr><tr><td>HIMGS 1100+2610</td><td>2021-03-09</td><td>7.6</td><td>1</td><td>low</td><td>23.4</td></tr><tr><td>HIMGS 1100+2621</td><td>2021-03-13</td><td>7.6</td><td>1</td><td>low</td><td>23.4</td></tr><tr><td>HIMGS 1100+2610</td><td>2021-03-14</td><td>7.6</td><td>1</td><td>low</td><td>23.4</td></tr></tbody></table>

**Note.** Column (1): The observation date. Column (2): The sample rate for recording the data. Column (3): The level of noise diode. Column (4): The rotated angle of the feed array.

slightly shorter than the integration time to avoid power leakage to the nearby time stamps. In addition, the noise diode spectrum is calibrated by observing the celestial point-source calibrator 3C286. The noise diode calibrator can be fired at either the high-power or the low-power level. With our test data, we found that the low-power level is sufficient for our calibration.

The observation of the 7 nights in 2021, that adopt the optimized observation parameter configuration, are the major survey data that will be used in the analyzes that follow. The pointing direction shifts by 10.835 per night in decl. direction. With the 7 nights drift-scan observations, a range of ~1° across the decl. direction is covered. With 4 hr drift scan, the observation covers the R.A. range from 9 to 13 hr, which overlaps with the Northern Galactic Cap (NGP) area of the Sloan Digital Sky Survey (SDSS; Reid et al. 2016). The detailed observation information is summarized in Table 1. The observation footprints of different data sets are shown in Figure 2.

### 3. Time-ordered Data Analysis

The raw data of the 19 feeds are dumped into files individually using the FITS<fnref n="15" /> (Flexible Image Transport System) format. Each FITS file contains a chunk of TOD with all frequency channels. The telescope-pointing direction data are recorded separately. To simplify the data analysis, we

convert the initial data format by combining the 19-feeds data as an extra axis and splitting the full frequency band into three sub-bands, i.e., the low-frequency band 1050–1150 MHz, mid-frequency band 1150–1250 MHz and high-frequency band 1250–1450 MHz. The telescope-pointing directions of the 19 feeds at each time stamp are calculated and written into the same data file.

The observation data include both the drift-scan data and the flux calibration data. The noise diode signal is used for real-time relative calibration of the gain, which is discussed in Sections 3.1 and 3.2. The noise diode flux spectrum is calibrated by observing a standard source, as discussed in Section 3.3. The RFI flagging is applied to the calibrated data. The details of the RFI flagging are described in Section 3.5. The data are then zero-centered by removing the temporal baseline variation before being used for mapmaking. The TOD data analysis pipeline is illustrated in Figure 3, where the blue rectangles indicate the input/output data and the green rectangles indicate the operation.

#### 3.1. Bandpass Calibration

The observed data value, $V(t, v)$, is the system gain multiplied with the combination of the input signal and noise

$$V(t, v) = g(t, v)(T(t, v) + n(t, v)), \quad (1)$$

where $T(t, v)$ is the antenna temperature corresponding to the total power collected by the telescope, $n(t, v)$ is the noise with $\langle n(t, v) \rangle_t = 0$, and $g(t, v)$ is the system gain. We assume that the

<footnote n="15">¹⁵ https://fits.gsfc.nasa.gov/</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 3. TOD analysis pipeline. The blue rectangles indicate the input/output data and the green rectangles indicate the operations performed on the data. The data preparation is described at the beginning of Section 3. The calibration pipeline includes bandpass calibration (Section 3.1), temporal drift calibration (Section 3.2), and flux calibration (Section 3.3). The main pipeline, the RFI flagging, baseline subtraction, and the mapmaking are described in Section 3.5, Section 3.4, and Section 3.6, respectively.</figcaption></figure>

system gain can be decomposed into a time-dependent component and a frequency-dependent component:

$$g(t, \nu) = g_t(t)g_\nu(\nu), \qquad (2)$$

where $g_\nu(\nu)$ is the bandpass gain factor and $g_t(t)$ is the temporal drift factor of the gain. In our analysis, $g_\nu(\nu)$ and $g_t(t)$ are calibrated using the noise diode, which is fired for 0.9 s in every 8 s as a relative flux calibrator. We assume that the noise diode temperature and spectrum are both stable during the observation. In every 8 s, we pick up the power value when the noise diode fired on, $V_{ND\ on}(t, \nu)$, and subtract the average power value at the two nearby time stamps, $V_{ND\ off}(t, \nu)$,

$$V_{ND\ on}(t, \nu) - V_{ND\ off}(t, \nu) = g_t(t)g_\nu(\nu)(T_{ND}(\nu) + n(t, \nu)), \qquad (3)$$

in which we assume both $g_t$ and the background emission are constant during the short time interval.

To check if the gain variation can be decomposed into factors of time variation and constant spectral shape as in Equation (2), we break the full drift scan into a few ~30 min time blocks ($t_b$) and evaluate the block-averaged bandpass gain

$$\bar{g}(t_b, \nu) = \frac{\langle V_{ND\ on}(t, \nu) - V_{ND\ off}(t, \nu) \rangle_{t_b}}{T_{ND}(\nu)}, \qquad (4)$$

where $\langle \cdots \rangle_{t_b}$ represents the averaging across the time block and $T_{ND}(\nu)$ is the noise diode spectrum. Because the noise diode was fired for ~1 s in every 8 s, the integration time within 30 min observation is 225 s. Assuming a system temperature of ~20 K, with 7.6 kHz frequency resolution, the system noise level is reduced to ~15 mK, which is negligible compared to the noise diode temperature of ~1 K.

The block-averaged bandpass is contaminated by the RFI. To remove the RFI contamination, we first smooth the data by

applying a median filter across the frequency channels to obtain an estimation of the bandpass baseline, and we then estimate the root mean square (rms) of the residual of the data after subtracting the smooth bandpass baseline. The frequency channels with values greater than 3 times the rms are flagged as RFI contaminated. We iterate the flagging until there are no extra masked channels. Note that the RFI flagging processing applied to the bandpass determination procedure is much stricter than that applied to the survey data. Some of the flagged channels here may actually not be real RFIs. However, it does not hurt to take a more strict criterion in the bandpass determination. Meanwhile, the median filter applied to the data is broad enough to recover a smooth bandpass baseline. After flagging all possible RFIs, the mask is then applied to the raw data without median-filter smoothing at each time stamp and the block-averaged bandpass is determined by taking the median value across the time block. A few of the frequency channels that were badly contaminated by RFI are fully flagged. The bandpass values at the fully flagged frequency channels are interpolated from the smoothed bandpass. Finally, to eliminate the noise, the bandpass is further smoothed with a third-order Butterworth low-pass filter with a critical delay frequency of $\tau = 0.7$ µs, which corresponds to a window size of 200 frequency bins. The choice of the window function size is further discussed in Section 4.1.

The overall gain drift of the block-averaged bandpass and the noise diode spectrum are removed by normalizing $\bar{g}(t_b, \nu)$ with its mean and producing the normalized bandpass

$$\bar{g}^n(t_b, \nu) = \bar{g}(t_b, \nu) / \langle \bar{g}(t_b, \nu) \rangle_\nu, \qquad (5)$$

where $\langle \cdots \rangle_\nu$ represents the averaging across the frequencies. To visualize the evolution of the shape of the bandpass, we check
<figure><img src="image_8.png" /><figcaption>**Figure 4.** The smoothed bandpass of the XX (left-hand panels) and YY (right-hand panels) polarizations for the 20210302 data, i.e., those taken on 2021 March 2. The raw frequency resolution is 7.6 kHz. Top panels: The block-averaged bandpass, each block is shown in a different color, and the black curve shows the averaged bandpass across all blocks. Bottom panels: The bandpass relative variation with respect to the first block. The gray areas mark the frequency channels with over 20% time stamps flagged for RFI.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** The relative variation of the bandpass between different nights with respect to the first night's observation (20210302).</figcaption></figure>

the normalized bandpass ratio with respect to the first block

$$ \gamma(t_b, \nu) = g^n(t_b, \nu) / g^n(t_0, \nu), \qquad (6) $$

where $t_0$ denotes the first time block.

An example of the results of Feed 01 for observation 20210302, i.e., those taken on 2021 March 2, is shown in Figure 4, and different time blocks are represented by different colors. The black curve is the mean bandpass across all time blocks. The gray areas mark the frequency channels with over 20% time stamps flagged. As shown by this figure, the bandpass has significant variation over frequency and also varies with time, but its shape is nearly constant over a few hours, except for a small fraction of frequency channels, which are badly contaminated by RFI. The bottom panels show the

bandpass ratio $\gamma(t_b, \nu)$, which has a variation less than 1% over 4 hr observation for most frequency channels.

Assuming that the shape of the bandpass is stable, Equation (3) is further averaged across the 4 hr of each night's observation

$$ \langle V_{\text{ND on}}(t, \nu) - V_{\text{ND off}}(t, \nu) \rangle_t = \bar{g}_t g_\nu(\nu) T_{\text{ND}}(\nu), \qquad (7) $$

where $\bar{g}_t = \langle g_t(t) \rangle_t$. The bandpass relative variation between the 7 nights observations in 2021 is shown in Figure 5. As an example, we show the two polarizations of the center feed in the left-hand and right-hand panels. The shape of the bandpass varies significantly between different days. Therefore, we emphasize that the shape of the bandpass needs to be determined for each night of observation.
<figure><img src="image_10.png" /><figcaption>**Figure 6.** The gain temporal variations of 4 hr observation. Each sub-panel shows the gain of one feed. The XX and YY polarizations are shown in red and blue colors, respectively. The dot markers represent the measurements with the noise diode; and the thick curves are the reconstructed gain variation.</figcaption></figure>

The bandpass-calibrated data is given by

$$V_t = \frac{V(t, \nu)}{\langle V_{\text{ND on}}(t, \nu) - V_{\text{ND off}}(t, \nu) \rangle_t} = \frac{g_t(t) T(t, \nu) + n(t, \nu)}{\bar{g}_t T_{\text{ND}}(\nu)} \quad (8)$$

## 3.2. Temporal Drift Calibration

The temporal drift is also calibrated with the noise diode. We average Equation (3) across frequencies

$$\langle V_{\text{ND on}}(t, \nu) - V_{\text{ND off}}(t, \nu) \rangle_{\nu} = \bar{g}_{\nu} g_t(t) (\langle T_{\text{ND}}(\nu) \rangle_{\nu} + n(t)), \quad (9)$$

where $\bar{g}_{\nu} = \langle g(\nu) \rangle_{\nu}$ and $n(t) = \langle n(t, \nu) \rangle_{\nu}$. If we further normalize with the time mean, then the temporal drift is

$$V_t = \frac{\langle V_{\text{ND on}}(t, \nu) - V_{\text{ND off}}(t, \nu) \rangle_{\nu}}{\langle V_{\text{ND on}}(t, \nu) - V_{\text{ND off}}(t, \nu) \rangle_{\nu, t}} = \frac{g_t(t)}{\bar{g}_t} + n'(t), \quad (10)$$

where we assume that both $\bar{g}_{\nu}$ and $T_{\text{ND}}$ are constant over time and $n'(t) = (g_t(t)/\bar{g}_t T_{\text{ND}})n(t)$. The first term of Equation (10), i.e., the normalized gain, represents the drifting of the actual gain, while the second term represents the variation caused by the measurement error in the calibration. $V_t$ represents the measurements of the gain value at each firing of the noise diode. Written in discrete form, we denote the gain measurements $V_t$ as vector $g_m$ and Equation (10) is expressed as

$$g_m = g_t + \mathbf{n} \quad (11)$$

We split the full-time stream into $\alpha$ short time blocks, $\Delta_{\alpha}$. The gain is assumed to be constant within each short block and varying between different blocks. Using a set of the base function $\mathbf{F} = \{F_1(t), F_2(t), \dots, F_{\alpha}(t)\}$, where

$$F_{\alpha}(t) = \begin{cases} 1 & t \in \Delta_{\alpha}, \\ 0 & \text{otherwise,} \end{cases} \quad (12)$$

the drifting of the gain is expressed as $g_t = \mathbf{F}g$, where $\mathbf{g}^T = \{g_1, g_2, \dots, g_{\alpha}\}$ is the parameter sets that need to be determined. In our analysis, we use a short block length of $\Delta_{\alpha} = 20$ s to avoid overfitting the temporal variation of the gain.

With the amplitude vector $\mathbf{g}$ as the parameter, and the measured gain values $g_m$, the likelihood is

$$P(\mathbf{g}_m) = P(\mathbf{g}_m | \mathbf{g}) P(\mathbf{g}) \propto \exp\left(-\frac{1}{2} \mathbf{n}^T \mathbf{N}^{-1} \mathbf{n}\right) \exp\left(-\frac{1}{2} \mathbf{g}^T \mathbf{C}_g^{-1} \mathbf{g}\right), \quad (13)$$

where $\mathbf{N}$ is the measurement noise covariance matrix, and $\mathbf{C}_g = \langle \mathbf{g} \mathbf{g}^T \rangle$ is the covariance matrix of the gain amplitude vector. The temporal variation can be modeled as follows:

$$\mathbf{C}_g = (\mathbf{F}^T \mathbf{F})^{-1} \mathbf{F}^T \mathbf{C}_N \mathbf{F} (\mathbf{F}^T \mathbf{F})^{-1}, \quad (14)$$

where the covariance matrix $\mathbf{C}_N$ is related to the noise power spectrum $P(f)$ as,

$$C_N(\delta t) = \int P(f) e^{2\pi i f \delta t} df. \quad (15)$$

Note that $P(f)$ only represents the power spectrum of the correlated noise (1/f noise), the total noise power spectrum is the combination of the white and correlated noise power spectrum, i.e., $P_{\text{total}}(f) = \frac{\sigma^2}{\delta\nu}(1 + P(f))$, where $\sigma$ is the systematic rms and $\delta\nu$ is the frequency resolution (Harper et al. 2018; Li et al. 2021).

The amplitude vector $\mathbf{g}$ can be solved by the maximum-likelihood method as

$$\hat{\mathbf{g}} = (\mathbf{F}^T \mathbf{N}^{-1} \mathbf{F} + \mathbf{C}_g^{-1})^{-1} \mathbf{F}^T \mathbf{N}^{-1} \mathbf{g}_m, \quad (16)$$

which is equivalent to Wiener filtering.

The measurements of $g_m$ for observation on 20210314, as an example, are shown in Figure 6 with the dot markers, the frequency range between 1150 to 1250 MHz is ignored due to the serious RFI contamination. The averaging bandwidth of two separated sub-bands is ~300 MHz in total. The measurements of different feeds are shown in different panels and the two polarizations are shown in red and blue colors, respectively. The corresponding temporal power spectrum is
<figure><img src="image_11.png" /><figcaption>Figure 7. The temporal power spectra of the time-ordered data. Each sub-panel shows the gain of one feed, the XX and YY polarizations are shown in red and blue colors. The dashed lines with markers show the measured power spectra and the solid lines show the best-fit model of Equation (17).</figcaption></figure>

shown with the dashed lines in Figure 7. There is a peak in the power spectrum at $f \sim 3 \times 10^{-4}$ Hz, corresponding to an oscillation in $V_t$ with period $\sim 1$ hour.<fnref n="16" /> Except for this, the power spectrum has the 1/f-type shape, which has higher power at the lower end of the f-axis. The 1/f-type shape power spectrum is due to the overall drift of the $V_t$ across time. The 1/f noise power spectrum is finally modeled as a combination of the 1/f-type power spectrum and a Lorenz profile,

$$ 
\begin{aligned}
P_{\text{total}}(f) &= \frac{\sigma^2}{\delta\nu}(1 + P(f)) \\
&= \frac{\sigma^2}{\delta\nu}\left(1 + \left(\frac{f}{f_k}\right)^\alpha + \frac{B}{1 + ((f - f_{\text{peak}})/w)^2}\right),
\end{aligned}
\quad (17) $$

where, $\alpha$, $f_k$, $B$, $f_{\text{peak}}$, and $w$ are the parameters that need to be fitted with the measured power spectrum. The solid lines in Figure 7 show the best-fit power spectrum.

With the best-fit noise power spectrum, we can estimate g and the temporal gain variation can be reconstructed with

$$ g_t = Fg. \quad (18) $$

The reconstructed temporal gains for observation 20210314 are shown with thick curves in Figure 6. The temporal gain variation is finally calibrated via

$$ V_2 = V_1/g_t. \quad (19) $$

Clearly, there are a few feeds, e.g., Feed 10 and Feed 17 of the data shown in Figure 6, that show significant gain fluctuation during the observation time. During the 7 nights of observations, these two feeds and Feed 06 sometimes show such gain fluctuation. This might be caused by hardware issues in the signal chains. This fluctuation increases the white noise level, i.e., the platform of the temporal power spectrum as shown in Figure 7, and is difficult to calibrate with our calibration strategy. We mark the feeds as “BAD FEED” when they are unstable and ignore them in the final mapmaking.

The temporal gain variations of the XX and YY polarizations are not congruous to each other, especially for Feed 08 and Feed 09 of the data shown in Figure 6. These large differences appear randomly in different feeds during the seven-night observations. Because the signal chains for the different polarizations are independent, we do not expect their gain variation to be identical or even similar to each other. As long as the drift of the gain is smooth and can be calibrated with our calibration strategy, the calibrated data are still used in the final mapmaking.

### 3.3. Absolute Flux Calibration

The absolute flux calibration is done by multiplying the temporal gain-calibrated data $V_2$ with the noise diode spectrum,

$$ T = V_2 \times T_{\text{ND}}, \quad (20) $$

where T is referred to below as the calibrated data, and $T_{\text{ND}}$ is referred to as the temperature of the noise diode. The noise diode temperature is measured via a series of hot-load measurements (Jiang et al. 2020) and it is assumed to be stable during the observations. During our observations, we performed several absolute flux calibrations using known celestial calibrators.

The absolute flux calibrations were made in drift-scan mode. The 19 feeds were grouped into five east-west lines. With five different pointings, the calibrator drifted across each feed in the same east-west line. To minimize systematic differences compared to the target field observation, a calibrator with its decl. close to the target field is required. We chose 3C286 as our flux calibrator and performed the calibration observation after the target observation of each day. The calibration pointing direction is shown in Figure 8.

The observation time is long enough to have the calibrator fully transit across the beam. The great-circle distance between the pointing direction and the calibrator, $\theta$, is calculated with

$$ \cos(\theta) = \sin(\delta_p)\sin(\delta_{\text{cal}}) + \cos(\delta_p)\cos(\delta_{\text{cal}})\cos(\alpha_p - \alpha_{\text{cal}}), \quad (21) $$

where $(\alpha_p, \delta_p)$ and $(\alpha_{\text{cal}}, \delta_{\text{cal}})$ are the R.A. and decl. of the pointing direction and the calibrator, respectively. We use the

<footnote n="16">¹⁶ The cause of the ~1 hr period oscillation is unknown. However, this oscillation is only observed in the 2021 data and it has disappeared in the later observations.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 8.** The footprint of the calibration observation. The red-dashed line represents 3C286’s drift route as it crosses the meridian in the local alt-az coordinate at FAST site. The calibrator drifts through each line of feeds with five different pointing directions. The circles packed in hexagons illustrate the position of the 19 feeds in each pointing direction.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** $\theta_{\text{FWHM}}$ measurements using 3C286 drift-scan observation. The circle markers show the mean of the best-fit $\theta_{\text{FWHM}}$ across the measurements of different days and the filled region indicates the corresponding rms. The results of different beam models are shown in different colors. The black squares show the measurements from Jiang et al. (2020).</figcaption></figure>

data within the time range with $\theta < 3'$ as the source-on power $V_{\text{on}}(\theta)$, and average across the time range with $\theta > 12'$ as source-off power $V_{\text{off}}$.

During the calibration, the noise diode is also fired in the same way as for the target field observations. The data are first calibrated against the noise diode power, $V_{\text{ND}}$, to cancel the bandpass gain. The corresponding main beam brightness temperature of the calibrator is

$$T_{\text{cal}} B(\theta) = \frac{V_{\text{on}}(\theta) - V_{\text{off}}}{V_{\text{ND}}} T_{\text{ND}}, \quad (22)$$

where $B(\theta) = \exp[-\frac{\theta^2}{2\sigma^2}]$ is the normalized beam pattern with $\sigma = \theta_{\text{FWHM}}/(2\sqrt{2\ln 2})$. The antenna temperature is converted from the source flux density via

$$T_{\text{cal}} = \frac{\eta_{\text{M}} \lambda^2}{2 k_{\text{B}} \Omega_{\text{MB}}} S, \quad (23)$$

where $k_{\text{B}}$ is the Boltzmann constant, $\lambda$ is the wavelength, $\Omega_{\text{MB}}$ is the main beam solid angle, and $\eta_{\text{M}}$ is the main beam efficiency. Assuming a symmetric Gaussian beam with half power beamwidth $\theta_{\text{FWHM}}$, the main beam solid angle is given by

$$\Omega_{\text{MB}} = 2\pi \left( \frac{\theta_{\text{FWHM}}}{2\sqrt{2\ln 2}} \right)^2 \approx 1.133\theta_{\text{FWHM}}^2. \quad (24)$$

The spectrum flux density of 3C286 can be modeled as in Perley & Butler (2017),

$$\log\left(\frac{S}{Jy}\right) = 1.2481 - 0.4507x - 0.1798x^2 + 0.0357x^3, \quad (25)$$

in which $x = \log(\frac{\nu}{\text{GHz}})$. The noise diode spectrum is evaluated by minimizing the following residual function for each feed, frequency, and polarization

$$\chi^2 = \sum_{\theta < 3'} \left| \frac{\lambda^2 B(\theta) S}{2 k_{\text{B}} \Omega_{\text{MB}}} - \frac{V_{\text{on}}(\theta) - V_{\text{off}}}{V_{\text{ND}}} \frac{T_{\text{ND}}}{\eta_{\text{M}}} \right|^2. \quad (26)$$

where $T_{\text{ND}}$ is the parameter to be determined.

We can also leave the half power beamwidth, $\theta_{\text{FWHM}}$, as another free parameter fit with the observation data. To model the sidelobes, we use the Jinc function beam model

$$B_{\text{Jinc}} = 4 \left( \frac{J_1(\pi\theta/\theta_{\text{FWHM}})}{\pi\theta/\theta_{\text{FWHM}}} \right)^2, \quad (27)$$

where $J_1(x)$ is the Bessel function of the first kind, and the cosine beam model

$$B_{\text{cos}} = \left( \frac{\cos(1.189\pi\theta/\theta_{\text{FWHM}})}{1 - 4(1.189\theta/\theta_{\text{FWHM}})^2} \right)^2, \quad (28)$$

which is known to have lower sidelobes compared to the Jinc function (Matshawule et al. 2021). We fit $\theta_{\text{FWHM}}$ at each
<figure><img src="image_14.png" /><figcaption>Figure 10. The best-fit η for different measurements.</figcaption></figure>

frequency with initial frequency resolution of 7.6 kHz. To reduce the variance, the best-fit values are then averaged in each 16 MHz frequency bin. The best-fits of the θ$_{FWHM}$ of the XX polarization from different beam models are shown in Figure 9 with different colors. The circle markers show the mean θ$_{FWHM}$ across the measurements in different days and the filled region indicates the corresponding rms. The beamwidth reported in Jiang et al. (2020) is shown with the black square markers. The best-fit θ$_{FWHM}$ of Feed 01, which is in the center of the FAST 19-feed array, shows consistent results across different days and beam models. However, the best-fit results of the other feeds reveal significant scattering between various days, e.g., Feed 02, or when using a different beam model, e.g., Feed 08 and Feed 18.

There are some tensions between our measurements of θ$_{FWHM}$ and the results of Jiang et al. (2020). For example, the θ$_{FWHM}$ of Feed 10 shows a significant deviation. Our results also show an obvious variation at the lower-frequency band. A couple of reasons could cause the differences in beam measurements. The beam shape might vary with the pointing direction because the reflector needs to be adjusted to form a paraboloid using different parts of the dish. During our drift-scan calibration observation, the telescope pointings are fixed in a few directions that are close to the Zenith of the FAST site. In this way, we could minimize the changes in the reflector. However, the drift-scan calibration is quite time-consuming and each beam can only be measured once during the transit of the celestial calibrator. With a single transit observation, we simply assume a symmetric beam profile and measure the beam profile across one direction for both the XX and YY polarizations. As shown in Jiang et al. (2020), the full beam shape is significantly asymmetric and can be well fitted using a “skew Gaussian” profile, which takes the ellipticity into account.

It is known that an accurate beam model is crucial for efficiently removing the foreground (Matshawule et al. 2021; Spinelli et al. 2022). With an accurate beam model, the systematic effects caused by the primary beam can be eliminated using the deep learning method (Ni et al. 2022; Gao et al. 2022). Current analysis indicates that the major primary beam effect of FAST is the asymmetric beam profile.

Symmetric beam models, such as the Jinc function or the cosine beam model, produce poor fitting results. The sidelobes of the FAST beam are poorly modeled using our current measurements. A complete analysis needs more observations and we will improve the measurements in further work. In the rest calibration and mapmaking analysis, we simply assume a Gaussian beam profile and adopt the beamwidth measurements of Jiang et al. (2020), and interpolate to the frequencies of our observation.

Because the noise diode signal is injected into the receiver system between the feed and low-noise-amplifier (LNA; Jiang et al. 2020), its sky-source-calibrated spectrum is slightly different from the noise diode spectrum model, which is measured using the hot-load. The difference is parameterized as

$$ \epsilon = \frac{T_{\text{ND}}^{\text{model}}}{T_{\text{ND}}}, \qquad (29) $$

where $T_{\text{ND}}^{\text{model}}$ is the noise diode spectrum model and $T_{\text{ND}}$ represents the noise diode spectrum determined using the celestial calibrator. In fact, $\epsilon$ degenerates with the main beam efficiency of FAST, $\eta_M$. Thus, $\epsilon$ and $\eta_M$ are combined as the total aperture efficiency $\eta = \eta_M\epsilon$ and determined using the celestial calibrator for each frequency, polarization, and feed.

It is known the aperture efficiency of FAST is weakly dependent on the Zenith Angle (ZA) with ZA < 26°4; and decreases quickly with ZA beyond 26°4 (Jiang et al. 2020). Because our calibration observations were always carried out when 3C286 is near its transit time, the pointing directions of the same feed are relatively consistent between different days and the ZA are all within 26°4. Thus, $\eta$ is assumed to be relatively stable between days for our calibration observations. The measured $\eta$ on different days are shown in Figure 10. Generally, the measured $\eta$ varies between different feeds but keeps a similar shape between different days for the same feed. The measured $\eta$ on each day is contaminated by RFI. To fill in the RFI gaps, we produce a template $\bar{\eta}$ for each feed by taking the median values of $\eta$ across the measurements on different days and fitting with a 15th-order polynomial function. The $\eta$ template is shown with the black-solid line in Figure 10. To recover the variations of $\eta$ between different days, the $\eta$
<figure><img src="image_15.png" /><figcaption>**Figure 11.** The temporal variation of baselines for the 6 days of observation. The different color shows different feeds. The temporal baselines for each feed are centered by subtracting the mean across the full observation, individually. The positive peaks are due to the bright continuum sources, and the few negative spikes are bad data due to the strong RFI contamination.</figcaption></figure>

measurement on each day is fitted to the template via

$$ \eta = (p_0 + p_1\nu + p_2\nu^2 + p_3\nu^3)\bar{\eta}, \quad (30) $$

where {$p_0, p_1, p_2, p_3$} are the parameters. The best-fit $\eta$ for different measurements is shown in Figure 10 with solid curves in the same colors as the corresponding measurements. The sky absolute flux density is finally obtained as

$$ T = V_2 \times T_{\text{ND}}^{\text{model}}/\eta. \quad (31) $$

### 3.4. Temporal Baseline Subtraction

We will call the average of the calibrated data across the frequency band the *baseline* of the data. The baselines for 6 days observation are shown in Figure 11 with different colors. The baselines are centered by subtracting the mean across the full observation time. The positive peaks in this data are due to bright continuum sources. Due to unknown reasons, some feeds occasionally perform badly during the observation, producing significantly larger fluctuations. Because different feeds point to different sky positions, we average the baselines across feeds to eliminate the flux variation from the sky. Such averaging across different feeds also reduces the baseline variance. Some of the data are contaminated by very strong RFI from satellites, which are fully flagged across the full frequency band. However, the data adjacent to these bad times are not flagged across the full frequency range. Since the baseline is estimated by taking the median value across the frequency band, these partially flagged time stamps may have significantly lower median values than the rest and result in negative spikes after subtracting the temporal mean.

In most cases, the baselines are steady, although the 20210306 data show some significant sharp variations, for which the reason is unknown. The shape of the baselines for different feeds is generally consistent, as expected, even for the 20210306 observation. This might indicate that the baseline

variation is due to a systematic background noise level variation during the observation time. We take the median across the 19 baselines to get rid of the spikes and further smoothed along the time with a median value filter. The smoothed mean baseline is shown with the thick black curve. The smoothed baseline is fitted to the calibrated TOD at each frequency

$$ A = (b^Tb)^{-1}b^T T, \quad (32) $$

where **b** represents the $n_t \times 1$ vector of baseline template, **T** represents the $n_t \times n_\nu$ matrix of TOD, and **A** is the $1 \times n_\nu$ vector of the fitting parameter. The baseline is subtracted via

$$ T_{\nu}^{c} = T_{\nu} - bA, \quad (33) $$

where $T^c$ represent the baseline centered TOD.

### 3.5. RFI Flagging

In Figure 12, the black curve shows the spectrum from one feed averaged over ~30 min. There is strong RFI contamination in the frequency band between 1150 to 1300 MHz, which are produced by the Global Navigation Satellite Systems (GNSSs), including the GPS, Galileo, Glonass, and Beidou (Teunissen & Montenbruck 2017). The frequency channels allocated for GNSS are marked in different colors. The contamination of such GNSS bands can leak into the neighboring channels due to the extended GNSS signal spectrum profile. We first remove the channels that are badly contaminated by GNSS, including those allocated frequency channels, as well as the neighboring channels within 1–20 MHz. The pre-removed frequency channels are shown with the gray area in Figure 12.

We then apply the SumThreshold and SIR (Scale-Invariant Rank) RFI flagging program (Offringa et al. 2010, 2012; Zuo et al. 2021) to the bandpass-calibrated data. The SumThreshold algorithm searches for consecutive
<figure><img src="image_16.png" /><figcaption>**Figure 12.** The spectrum from a feed averaged over ~30 min. The channels contaminated by RFI from GNSS are marked by the gray region. The allocated frequency channels used for the GPS, Galileo (GAL), Glonass (GLO), and Beidou (BD2/BD3) are marked in different colors (Teunissen & Montenbruck 2017).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Waterfall plot of ~4 hr TOD observed with Feed01 XX polarization. The blank region is the RFI-flagged data.</figcaption></figure>

points of different numbers (increase from 1 to $2^n$) in TOD as potential RFI contaminated points that have a value above certain preset thresholds, i.e., $\chi_n = \chi_1/1.5^{\log_2 n}$, where $\chi_1 = 10$ is the initial threshold and $n$ is the number of data sample considered. The thresholds vary according to the number of data samples considered. To find the extra weak contamination near the flagged high values, the SIR RFI flagging is then applied. The SIR RFI flagging method uses the one-dimensional mathematical morphology technique to find the neighbors intervals in the time or frequency domain that are likely to be affected by RFI.

Compared to other radio telescopes, FAST is much more sensitive, with many genuine celestial radio sources, and even the HI emissions from nearby galaxies can be detected with a high signal-to-noise ratio (S/N) in the raw data; therefore, care must be taken to avoid removing them by mistake. We stack the same time outputs of the different feeds, which would lower the celestial source signal because the different feeds are pointed at slightly different sky directions at any given time while enhancing the RFIs that enter through the beam side lobe and are simultaneous on all feeds. In addition, the Sum-Threshold algorithm is only applied along the frequency axis.

As an example, the RFI flagged 20210314 data for Feed01 XX polarization is shown in a waterfall plot in Figure 13. The blank regions are flagged data. Most frequency bands between 1150 and 1323 MHz are flagged because they are badly contaminated by the GNSS. Another severely RFI

contaminated part is around 1090 MHz, which is the band allocated to aircraft Automatic Dependent Surveillance-Broadcast (ADS-B) data communication, where the RFI occurs frequently during the 4 hr observation. The frequency range beyond 1323 MHz is relatively free of RFI.

To check the RFI residual after flagging, we show a histogram of the data. The contribution of the natural continuum emission can be removed by taking the difference between the two neighboring frequency channels. The histograms are produced using data in three different frequency ranges, i.e., 1050–1150 MHz, 1150–1323 MHz, and 1323–1450 MHz, and the results are shown in Figure 14. The results of the data before RFI flagging are shown with solid lines and those after RFI flagging are shown with dashed lines. Before RFI flagging, the data shows a combined profile of a Gaussian distribution and a high-temperature tail, which indicates a significant RFI contamination. After the RFI flagging, the high-temperature tails are greatly reduced in all three frequency bands. The histogram statistic also shows that, with the SumThreshold flagging, there are 34.7%, 79.2%, and 12.2% data flagged in the three frequency bands, respectively; and with the additional SIR flagging, another 6.5%, 8.5%, and 2.0% data are flagged. For the data on other days, the ratios are more or less similar.

The frequency channels within 1150–1323 MHz are badly contaminated by the strong RFI contamination. In particular, the channels at the lower end of the frequency band, 1250–1450 MHz, contain some GNSS signal bands. Therefore,
<figure><img src="image_18.png" /><figcaption>Figure 14. Histogram statistic of the temperature difference before and after the RFI flagging.</figcaption></figure>

we adjust the valid frequency range of the high-frequency band to 1323–1450 MHz.

### 3.6. Maps

The calibrated data are zero-centered by subtracting the baseline, and the XX and YY polarization are combined into Stokes I. The TOD is then projected to the map domain via the standard mapmaking procedure (Tegmark 1997). To save the computation time for mapmaking, we re-binned the data to 28 kHz frequency resolution. The maps are made for each frequency without considering the correlation between different frequency channels.

We use the variance across time of the ~30 min data block for each feed and polarization as the noise variance. The noise is assumed to be uncorrelated between different time blocks and its covariance matrix is assumed to be diagonal, i.e., $N = \text{diag}\{\sigma^2(t)\}$. The map is obtained by

$$ \hat{\mathbf{m}} = (\mathbf{P}^\mathrm{T}\mathbf{N}^{-1}\mathbf{P})^{-1}\mathbf{P}^\mathrm{T}\mathbf{N}^{-1}\mathbf{T}, \quad (34) $$

in which $\mathbf{P}$ is the pointing matrix, which relates the time to the map coordinate. We use the HEALPix scheme for the sky with NSIDE = 2048, corresponding to a pixel size of 1'72 (pixel area of 2.95 arcmin²).

Figure 15 shows the frequency-averaged maps. The maps made with the data of two RFI-free frequency bands are shown in the top and bottom panels, respectively. Since the surveyed sky is a long strip, we divide the full region into several pieces in the R.A. direction, each is about 1 hr (i.e., 15°) in R.A. We can see many point sources on the map clearly, and the point sources in the 1050–1150 MHz band have good correspondence with the point sources in the 1323–1450 MHz band, which is what we would expect for radio continuum sources such as quasars and radio galaxies.

To improve the flux measurements of point sources, we applied an alternative mapmaking procedure similar to Haynes et al. (2018)

$$ \hat{\mathbf{n}}_p = [(\mathbf{P}\mathbf{K})^\mathrm{T}\mathbf{T}]_p / [(\mathbf{P}\mathbf{K})^\mathrm{T}\mathbf{I}]_p, \quad (35) $$

where $p$ indicates the $p$th pixels of the map, $\mathbf{I} = \{1, 1, \dots, 1\}^\mathrm{T}$ is a column vector with all elements equal to 1, and $\mathbf{K}$ represents a kernel function that paints the antenna temperature to the

nearby pixels. We use a Gaussian kernel function

$$ K_{pq} = \exp \left[ -\frac{1}{2} \left( \frac{r_{pq}}{\sigma_K} \right)^2 \right], \quad (36) $$

where $K_{pq}$ is the element of the $\mathbf{K}$ matrix, $r_{pq}$ is the great-circle distance between the $p$th and $q$th pixels of the map, and $\sigma_K$ indicates the kernel size. We set $\sigma_K = 1.5$ and use HEALPix scheme with NSIDE = 4096, corresponding to a pixel size of 0'86. To make the mapmaking process easier, the TOD from 1375 to 1425 MHz are averaged into a single frequency channel before the mapmaking. This map is used for flux comparison with the continuum sources. We extract the flux of point sources and compare them with the NVSS continuum measurements. Finally, 81 isolated point sources with flux over 14 mJy are identified within the surveyed region. The flux measurements are presented in Table 3 and the detailed discussion can be found in Section 4.5.2.

Since the observed data is already the convolution of the sky signal and the telescope beam pattern, using the kernel function Equation (36) to create the map is equivalent to an additional convolution. If we assume a Gaussian beam model with a beamwidth of $\theta_{\mathrm{FWHM}}$, the final map is then smoothed with a Gaussian function with a kernel size of $\sigma_{K'} = \sigma_K + \theta_{\mathrm{FWHM}}/(2\sqrt{2\ln 2})$.

## 4. Discussion

### 4.1. Bandpass Ripple

The bandpass calibration changes the shape of the spectrum significantly. To understand the ripple structure in the bandpass, we estimate the delay spectrum $|\tilde{V}(\tau)|^2 = \langle \delta\tilde{V}(\tau)\delta\tilde{V}^*(\tau) \rangle$, where

$$ \delta\tilde{V}(\tau) = \int d\nu \, \delta V(\nu) e^{-2\pi i\nu\tau} $$

is the Fourier transform of the data across the frequency. We make the analysis with the data contrast, i.e., $\delta V(\nu) = \frac{V(\nu)}{\bar{V}} - 1$, where $V(\nu) = \langle V(t, \nu) \rangle_t$ is the data value at frequency $\nu$ averaged across observation time and $\bar{V} = \langle V(\nu) \rangle_\nu$ is the mean $V(\nu)$ averaged across the frequency band. The delay spectrum is taken for the data before the bandpass calibration, after bandpass calibration, and also for the measured bandpass itself. The bandpass ripple structure would show up as a peak in the delay spectrum at a particular delay value.

In Figure 16, we plot the delay spectra of the bandpass (top panel) and the data (bottom panel) of the XX polarization of the center feed during the 20210302 observation. The spectra for the other feeds and polarizations are similar. For the bandpass (top panel), as mentioned in Section 3.1, we show the results of smoothing with the third-order Butterworth low-pass filter of three different sizes (defined as the 3 dB compression point) $\tau_c = 1.3$ µs, 0.7 µs, and 0.4 µs, respectively. As expected, the delay spectra are strongly suppressed at the small scale by the smaller-sized filters. For the delay spectra, we show the raw data (i.e., TOD before calibration) and the calibrated data using the bandpass smoothed with the three different filter sizes.

We can see an obvious peak in the bandpass delay spectrum, which is marked by a dashed vertical line in the figure. This peak is associated with a standing wave between the FAST
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Total intensity map of the surveyed regions. The top and bottom panels show the integrated intensity between 1050 and 1150 MHz and 1323–1450 MHz, respectively.</figcaption></figure>

feed and reflector, with a delay of $\tau \sim 0.92$ µs, corresponding to a standing wave with a ripple wavelength of 1.087 MHz in the spectrum. However, we do not see a significant standing-wave peak in the raw TOD, the standing-wave signature is more prominent in the bandpass measurements than in the sky observation, probably because this standing wave is induced by the noise diode itself. If the data is calibrated with such a bandpass, it would induce the ripple structure that is not present in the raw data itself. This can be avoided by employing the

bandpass smoothed with small-sized filters, e.g., those with $\tau_c = 0.7$ µs and $\tau_c = 0.4$ µs filters, as shown in the bottom panel of Figure 16. Note that the smoothing filter is applied to the bandpass and not the sky spectrum.

The delay spectra of the TOD before bandpass calibration for all 7 nights are shown in the top panel of Figure 17. Two out of the 7 nights' data show weak standing-wave peaks. The standing-wave signature is generally consistent across different polarization and beams within the same night's observation.
<figure><img src="image_20.png" /><figcaption>**Figure 16.** Top panel: The delay spectrum of the bandpass. The delay spectra of the bandpass smoothed with different window function sizes are shown in different colors. Bottom panel: The delay spectrum of the TOD. The gray curve shows the delay spectrum of the TOD before bandpass calibration, while the other curves show the results of bandpass-calibrated data using bandpass measurement smoothed with different window function sizes. All of the delay spectra are estimated using the XX polarization of the central beam data of 20210302. The dashed vertical line marks the delay of $\tau \sim 0.92$ µs, associated with the standing wave between the FAST feed and reflector.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 17.** The TOD delay spectrum of different days. The top panel shows the delay spectra of the TOD before bandpass calibration, while the bottom panel shows the results after bandpass calibration. All the delay spectra are estimated using the XX polarization of the central beam data. The dashed vertical line marks the delay of $\tau \sim 0.92$ µs, associated with the standing wave between the FAST feed and reflector.</figcaption></figure>

We use a low-pass filter with a size of $\tau_c = 0.7$ µs to suppress both the noise and standing-wave signature during the bandpass determination. The delay spectra of the bandpass-calibrated data are shown in the bottom panel of Figure 17, the standing-wave peak is negligible.

<figure><img src="image_22.png" /><figcaption>**Figure 18.** The histogram of the residual TOD. The solid green curve represents the combined data histogram and the dashed-red line represents the best-fit Gaussian distribution. The gray-thin curves represent the histograms of different data blocks. Each curve is normalized by its total number of data points.</figcaption></figure>

## 4.2. Measurement Uncertainty

The system temperature, $T_{sys}$, is related to the measurement noise level via the radiometer equation

$$ \sigma = \frac{T_{sys}}{\sqrt{N_{pol}N_{feed}\Delta t\Delta\nu}}, \qquad (37) $$

where $\sigma$ is the rms of the measurements representing the noise level, $N_{pol} = 2$ is the number of polarization, $\Delta t = 1$ s is the integration time, $\Delta\nu = 28$ kHz is the frequency resolution, and $N_{feed}$ is the number of the feed. We estimate the measurement noise level for each of the feeds individually and adopt $N_{feed} = 1$ in the following analysis.

To calculate the rms, we subtract the continuum emission of the point sources. We revise the rms estimation method introduced in Wang et al. (2021), which uses the difference between four adjacent frequency channels

$$ \Delta T(\tilde{\nu}) = \frac{1}{2}(T(\nu_1) + T(\nu_3)) - \frac{1}{2}(T(\nu_2) + T(\nu_4)), \quad (38) $$

where $\nu_1 < \nu_2 < \nu_3 < \nu_4$ are the four adjacent frequency channels and $\tilde{\nu} = (\nu_1 + \nu_2 + \nu_3 + \nu_4)/4$ is the reduced frequency of the residual data. This can remove most of the continuum emission from the sky that is linear across the frequencies. The rms of the residual data is related to the original data rms via

$$ \sigma(\tilde{\nu}) = \sqrt{\frac{1}{4}(\sigma^2(\nu_1) + \sigma^2(\nu_3)) + \frac{1}{4}(\sigma^2(\nu_2) + \sigma^2(\nu_4))}. \quad (39) $$

We calculate the noise level using the calibrated data. Each day's data is divided into seven blocks, each ~30 minute. The two polarizations are combined into the total intensity, i.e., $T = (T_{XX} + T_{YY})/2$. Figure 18 shows the histogram statistics of the residual data value. For each time block, the histogram statistic includes all of the 19 beams' data and the result is shown with the gray curve. The histogram for all of the time blocks combined is shown with the green curve. All of the histograms are normalized with the total number of data samples. The results of the low-frequency band, i.e., 1050–1150 MHz, are shown in the left-hand panel, and the
<figure><img src="image_23.png" /><figcaption>Figure 19. The noise level of each feed. The results for lower and higher frequency bands are shown in red and blue markers. The feeds are split into three categories, i.e., the feeds in the central (Feed 01), the inner circle (Feed 02–Feed 07), and the outer circle (Feed 08–Feed 19) of the feed array. The horizontal lines indicate the mean noise level of each feed category.</figcaption></figure>

results of the high-frequency band, i.e., 1323–1450 MHz, are shown in the right-hand panel, respectively.

We fit the histogram with a Gaussian function,

$$N(T) = A \exp \left[ -\frac{1}{2} \frac{(T - \mu_0)^2}{\sigma^2} \right]. \quad (40)$$

The best-fit function of the all-time-block combined histogram is shown with the dashed-red curve in Figure 18 and the best-fit $\sigma$ and $\mu_0$ are also shown in the legend. Both of the two frequency bands’ data fit the Gaussian function well, which indicates that the residual data are dominated by white noise. The best-fit noise levels of the data are 116.39 mK and 99.47 mK for the low-frequency band and high-frequency band, respectively.

The system temperature includes several different components and can be expressed as,

$$T_{\text{sys}} = T_{\text{rec}} + T_{\text{sky}} + T_{\text{CMB}}, \quad (41)$$

where $T_{\text{rec}}$ is the receiver temperature, $T_{\text{sky}}$ is the sky temperature, and $T_{\text{CMB}} = 2.725$ K is the mean brightness temperature of the cosmic microwave background (CMB). We ignore the temperature from ground-spill when the telescope is pointing close to the Zenith. The major component of the $T_{\text{sky}}$ is the diffuse emission of the Galactic synchrotron. According to the measurements in the work of Jiang et al. (2020), the system temperature close to the Zenith is about 20 K. Using Equation (37) and substituting $N_{\text{pol}} = 2$, $N_{\text{feed}} = 1$, $\Delta t = 1$ s, and $\Delta \nu = 28$ kHz, we should have $\sigma = 84.5$ mK. The slightly higher noise level could be caused by residual continuum emissions.

The noise level also fluctuates between different feeds. As shown in Figure 19, the noise level of the lower and higher frequency bands are shown in red and blue markers. The feeds are split into three categories, i.e., the feeds in the central (Feed 01), the inner circle (Feed 02–Feed 07), and the outer circle (Feed 08–Feed 19) of the feed array. The horizontal lines indicate the mean noise level of each feed category. Although the noise levels of different feeds are varying, there is a trend for the feeds in the outer circle of the feed array to have higher noise, i.e., about 10% increases in the noise level with respect to the central feed.

This system temperature estimation can also be applied in the map domain. We estimate the pixel noise level of the map at each reduced frequency, $\tilde{\nu}$. The histogram of each reduced

<figure><img src="image_24.png" /><figcaption>Figure 20. Same as Figure 18 but for the residual maps. The gray-thin curves represent the histograms of different frequencies.</figcaption></figure>

frequency is shown with the gray curve in Figure 20. The green curve shows the total histogram using all frequencies. The results of the low/high-frequency band are shown in the left-hand and right-hand panels. The dashed-red curve shows the best-fit Gaussian function. Only histograms with amplitudes greater than half of their maximum are used in the Gaussian function fitting. The best-fit $\sigma$ of the Gaussian function indicates the pixel noise level.

Clearly, the total histogram profile departs from the Gaussian function at $|T| \gtrsim 0.05$ K. The pixel noise levels for the low-frequency and high-frequency bands are 36.0 mK and 25.9 mK, respectively. Because the fit only uses the histograms with amplitudes greater than half of their maximum, the pixel noise level does not take into account the effect of the large residual values. Nevertheless, the pixel noise level is higher than the forecast. Due to the different RFI flagging fractions, the mean integration times for the low-frequency and the high-frequency band are 10.0 s and 16.7 s, respectively. Thus, assuming $T_{\text{sys}} = 20$ K, the pixel noise levels for these two frequency bands are 26.8 mK and 20.7 mK, respectively. Several reasons could potentially increase the noise; for example, the weak RFI contamination, which is below the noise level of the original TOD, becomes dominant when the pixel noise level is lowered by integrating data via the mapmaking process; the residual sky contamination is not removed with Equation (38); or there is weakly correlated noise in the original TOD.

### 4.3. RFI Residual

As discussed in Section 3.5, most of the RFI contamination comes from the signal transmission of GNSS. Even after the RFI flagging, there may be some RFI residuals. For instance, as shown in Figure 14, there are still tails at the large-value-end deviating from the Gaussian distribution. The non-Gaussian tail of Figure 20 may also be due to the RFI residuals.

These residuals are mostly weak and probably exist in the frequencies close to the flagged channels. A more strict flagging may reduce the residuals. However, it might also cause significant data loss. In particular, it may cause over-flagging of the HI emission line of galaxies at lower redshift. Although this misidentification can be partly avoided by stacking the same time outputs of the different feeds, some galaxies may be bright enough to be flagged as RFI.

Instead of performing a stricter RFI flagging to the TOD, it is better to mitigate the residuals after making the map. For instance, in the analysis of Cunnington et al. (2023), a final aggressive RFI
<table><caption>**Table 2**

Bright Sources Used for Bandpass Shape Inspection</caption><thead><tr><th>Source Name</th><th>74 MHz<sup>a</sup></th><th>151 MHz<sup>b</sup></th><th>365 MHz<sup>c</sup></th><th>408 MHz<sup>d</sup></th><th>1400 MHz<sup>e</sup></th><th>4850 MHz<sup>f</sup></th></tr></thead><tbody><tr><td>B2 1010+27A</td><td>1440 ± 170</td><td>910 ± 66</td><td>424 ± 26</td><td>470 ± 72</td><td>158 ± 4.8</td><td>40 ± 6.0</td></tr><tr><td>B2 1016+26</td><td>...</td><td>700 ± 63</td><td>357 ± 25</td><td>304 ± 70</td><td>121 ± 4.3</td><td>39 ± 5.9</td></tr><tr><td>B2 1018+26</td><td>1660 ± 190</td><td>1400 ± 79</td><td>640 ± 57</td><td>735 ± 80</td><td>228 ± 7.7</td><td>76 ± 11</td></tr><tr><td>B2 1036+26</td><td>1180 ± 140</td><td>870 ± 59</td><td>...</td><td>354 ± 70</td><td>145 ± 4.7</td><td>46 ± 6.9</td></tr><tr><td>B2 1037+26</td><td>1840 ± 260</td><td>860 ± 63</td><td>408 ± 32</td><td>425 ± 72</td><td>124 ± 4.4</td><td>30 ± 7.0</td></tr><tr><td>B2 1039+27</td><td>1880 ± 200</td><td>1200 ± 66</td><td>559 ± 30</td><td>603 ± 72</td><td>149 ± 4.5</td><td>31 ± 4.7</td></tr><tr><td>B2 1054+26</td><td>4620 ± 490</td><td>2650 ± 133</td><td>1520 ± 53</td><td>1260 ± 100</td><td>436 ± 15</td><td>95 ± 14</td></tr><tr><td>B2 1138+27</td><td>2150 ± 230</td><td>1380 ± 97</td><td>816 ± 26</td><td>873 ± 80</td><td>284 ± 8.5</td><td>92 ± 14</td></tr><tr><td>B2 1202+26</td><td>2250 ± 240</td><td>2110 ± 124</td><td>1250 ± 172</td><td>969 ± 80</td><td>379 ± 11</td><td>141 ± 21</td></tr><tr><td>B2 1203+26B</td><td>770 ± 130</td><td>670 ± 98</td><td>415 ± 38</td><td>420 ± 72</td><td>156 ± 4.7</td><td>63 ± 9.5</td></tr></tbody></table>

**Notes.** The flux values are all in the unit of mJy.

<sup>a</sup> The VLA Low-Frequency Sky Survey (Cohen et al. 2007).

<sup>b</sup> The 7C survey of radio sources at 151 MHz (Waldram et al. 1996).

<sup>c</sup> The Texas Survey of Radio Sources (Douglas et al. 1996).

<sup>d</sup> The B2 Catalogue of radio sources (Colla et al. 1972).

<sup>e</sup> The NRAO VLA Sky Survey (Condon et al. 1998).

<sup>f</sup> A New Catalog of 53522 4.85 GHz Sources (Becker et al. 1991).

<figure><img src="image_25.png" /><figcaption>**Figure 21.** The spectra of bright sources in the field. The red markers show the measured source spectrum. The flux is averaged within each 10 MHz frequency bin and the error bar indicates the flux rms in each frequency bin. The gap between 1150 and 1300 MHz is due to RFI contamination. The dashed line indicates the source spectrum model fit using flux measurements in the literature and the gray area indicates the spectrum model uncertainty.</figcaption></figure>

flagging is performed on the eigenmodes of the principal component analysis. For the case of HI galaxy searching, the RFI residuals can be identified by comparing the measurements between different feeds because the RFI residuals are supposed to be shown in more than one feed simultaneously.

Additionally, there might be another kind of RFI residual, as pointed out by the GNSS simulation work of Harper & Dickinson (2018). The GNSS signal may have out-of-band leakage to a broad frequency range. The leakage signal of all of the satellites close to the point direction contributes an extra noise, which may be comparable to the equipment noise level. The noise level of the leakage is also affected by the beam's side lobe level and can be evaluated with a detailed simulation in future work.

### 4.4. Spectra of Sources

We also inspect the bandpass shape by comparing the spectra of bright sources to the flux measurements in the

literature. We chose the bright sources that are closely scanned by at least one beam, i.e., the minimal angular distance between the source and the beam center is less than 0.5. Meanwhile, the source should have flux measurements at multiple frequency bands. In this analysis, we use 10 radio sources with flux measurements at 74 MHz (Cohen et al. 2007), 151 MHz (Waldram et al. 1996), 365 MHz (Douglas et al. 1996), 408 MHz (Colla et al. 1972), 1.4 GHz (Condon et al. 1998), and 4.85 GHz (Becker et al. 1991). The flux of these 10 sources is listed in Table 2. The source spectrum is modeled by fitting a third-order polynomial function to the flux measurements.

The source spectra are extracted from the calibrated data by taking the spectra at the time stamp when the source center is mostly close to the pointing direction. The data within the frequency band 1150–1250 MHz is ignored due to serious RFI contamination. The spectrum of each source is then averaged in every 10 MHz frequency bin. The measured spectra of the 10 sources are shown in Figure 21. The error bar indicates the rms
<figure><img src="image_26.png" /><figcaption>**Figure 22.** The differential number count of the continuum sources in the survey area. The black stepping curves display the differential number count of the NVSS catalog, while the red circle markers with solid lines display the results of our map. The green vertical dashed-dotted line indicates the flux limit of $5\sigma_c \approx 7$ mJy beam$^{-1}$, where $\sigma_c$, shown with the green vertical dashed line, is the confusion limit due to the background unresolved sources. The gray vertical dashed line indicates the 30 $\mu$Jy beam$^{-1}$ flux limit of the map.</figcaption></figure>

of the flux measurements within each 10 MHz frequency bin. The polynomial-fitted source spectrum model is shown with the black-dashed line and the gray area indicates the model uncertainty, i.e., the upper/lower bound is estimated by fitting the third-order polynomial function to the upper/lower limit of 68% flux measurement confidence interval. The gap between 1150 and 1300 MHz is due to RFI contamination.

Generally, our measurements produce a smooth power-law shape spectrum, which indicates that the bandpass calibration efficiently corrects the bandpass shape. The spectrum shape slightly fluctuated at frequencies close to the RFI contamination, especially for the relatively faint sources. The flux is generally consistent with the spectrum model fitted using the flux measurements at a few frequency bands in the literature. The deviation between our measurement and the spectrum model, e.g., source B2 1039+27, might happen because the intrinsic spectrum variation of the source cannot be well-fitted by a low-order polynomial function.

## 4.5. Flux of Detected Sources

According to the measurements in Section 4.2, in our survey the pixel noise level is ~25.9 mK with a frequency resolution of 28 kHz at high-frequency band. The corresponding flux limit at 1400 MHz with a bandwidth of 50 MHz should be ~30 $\mu$Jy beam$^{-1}$. We apply a source-finding algorithm, i.e., the *DAOStarFinder*,<fnref n="17" /> to our map with the threshold of 30 $\mu$Jy beam$^{-1}$ and aperture size of 3' and more than three thousand continuum sources are detected. The flux-weighted differential number count of the detected continuum sources is shown using the red circle markers with the solid curve in Figure 22.

However, we should also consider the confusion limit for the continuum sources (Condon 1974; Meyers et al. 2017),

$$ \sigma_c \approx 0.2 \left( \frac{\nu}{\text{GHz}} \right)^{-0.7} \left( \frac{\theta_{\text{FWHM}}}{\text{arcmin}} \right)^2 \approx 1.4 \text{ mJy beam}^{-1}. \quad (42) $$

This is much larger than the limit given earlier. Sources fainter than a few $\sigma_c$ would be confused and not detected as individual sources.

To check our survey results, we compare the continuum flux density of the sources in the observed field with those in the NRAO-VLA Sky Survey (NVSS) catalog (Kimball & Ivezić 2008). We use the integrated flux density of NVSS sources from a combined radio objects catalog with flux and position corrections.<fnref n="18" /> The flux limit of the NVSS catalog is given as 2 mJy beam$^{-1}$. There are 3161 NVSS sources in the survey area, i.e., $9^h < R.A.< 13^h$ and $+25^\circ 83 < Dec < +27^\circ 08$. The flux-weighted differential number count for sources in the NVSS catalog is also shown in Figure 22 with the black stepping curves. The detected continuum sources using our map is consistent with the NVSS down to ~7 mJy beam$^{-1} \approx 5\sigma_c$. At the faint end below 7 mJy (marked in the figure by the vertical dashed-dotted line), the number of sources detected by our survey begins to fall below that of the NVSS, which is unsurprising because the NVSS has much higher angular resolution, and therefore a lower confusion limit.

To make source-by-source flux measurement comparison, we select isolated bright sources from the full NVSS sample according to the following criteria:

1. We reject the sources that have neighbors' flux over 10% of the central source within $9'$, i.e., about three times of the beamwidth (Gregory et al. 1996). This selection criteria rejects more than 90% of the NVSS sources in the field.

2. We then remove the sources with flux less than 14 mJy beam$^{-1}$. We adopt such an aggressive flux limit to avoid confusion from the background noise. Another 3% source is rejected according to this criteria.

3. Finally, we pick the closely scanned sources that are $0.5\theta_{\text{FWHM}} \sim 1.5$ or less from the center of at least one FAST beam.

A total of 81 isolated sources meet these selection criteria, making up the isolated sample. This sample is listed in Table 3, which is given in the appendix. This isolated sample of sources is used for the following source-by-source flux measurement comparison.

### 4.5.1. Flux Comparison with Time-ordered Data

We first make a comparison of flux from the TOD. We use the mean flux density across the frequency range 1.375–1.425 GHz, which is the same frequency range of the NVSS catalog (Condon et al. 1998). Because the sky coverage partially overlaps between different days, the same source may be observed by different beams on different days. The number of sources used for each feed is listed below,

<table><tr><td>Feed #</td><td>01</td><td>02</td><td>03</td><td>04</td><td>05</td><td>06</td><td>07</td><td>08</td><td>09</td><td>10</td></tr><tr><td>N</td><td>19</td><td>24</td><td>17</td><td>23</td><td>9</td><td>15</td><td>18</td><td>18</td><td>12</td><td>5</td></tr><tr><td>Feed #</td><td>11</td><td>12</td><td>13</td><td>14</td><td>15</td><td>16</td><td>17</td><td>18</td><td>19</td><td></td></tr><tr><td>N</td><td>19</td><td>8</td><td>20</td><td>19</td><td>23</td><td>15</td><td>10</td><td>23</td><td>10</td><td></td></tr></table>

The measured flux density is extracted by taking the spectrum density at the time when the source has the minimal angular distance to the feed center and is compared with the expected flux, which is obtained by multiplying The NVSS flux density with a Gaussian beam profile according to the angular

<footnote n="17">¹⁷ https://photutils.readthedocs.io/en/stable/api/photutils.detection.DAOStarFinder.html</footnote>

<footnote n="18">¹⁸ http://www.aoc.nrao.edu/~akimball/radiocat.shtml</footnote>
distance to the beam

$$ \hat{S}_{\text{NVSS}} = S_{\text{NVSS}} \exp \left[ -\frac{1}{2} \frac{\theta_{\min}^2}{(\theta_{\text{FHWM}}/(2\sqrt{2\ln 2}))^2} \right]. \quad (43) $$

The flux–flux comparison is plotted in Figure 23. The NVSS sources scanned by different feeds are shown with different colors and those sources scanned by the central, inner circle, and outer circle of the FAST feed arrays are shown in the top, middle, and bottom panels, respectively. The measured flux densities are shown to be consistent with the NVSS catalog.

We quantified the scattering of the measurements using the relative flux residual with respect to the total flux

$$ \delta S = (S - \hat{S}_{\text{NVSS}}) / \sqrt{S \hat{S}_{\text{NVSS}}}, \quad (44) $$

where S represents the extracted fluxes from our TOD. The histogram statistic of the flux residuals is illustrated in Figure 24. The sources scanned by different feeds are also shown with different colors and those sources scanned with the feeds in the central, the inner circle and the outer circle of the feed array is shown in the top, middle, and bottom panels, respectively. The black-solid curves show the averaged histogram of the feed categories. All of the histograms are normalized with the total number of measurements. The rms of the relative flux residual, i.e., $\sqrt{\sum \delta S^2 / N_{\text{total}}} \times 100\%$, of the three feed categories are 4.2%, 5.3%, and 6.9%, respectively. The source flux measured with the central feed has less scattering than measurements with the other feeds. The increasing residue for feeds in the inner and outer circle of the feed array is probably a result of an error in the beam model because the beams are more distorted as we move out from the center.

We also check the flux measurements uncertainty between the observation on different days. The results show that the flux residual rms of different days are generally consistent. We can also take advantage of the repeated observation of the same strip on 2021 March 9, and 2021 March 14. We estimate the relative residual rms using the flux differences of the same sources between these two observations.

$$ \sqrt{\sum \delta S_{\text{ab}}^2 / N_{\text{total}} / 2} \times 100\% = 3.7\%, $$

where $\delta S_{\text{ab}} = (S_a - S_b) / \sqrt{S_a S_b}$ is the flux difference. Since the observations on these two days have the same pointing direction, the systematic effect, such as the beam effect, is canceled. If the source flux variation between the short period is negligible, then this relative residual rms indicates the calibration uncertainties in our point-source flux measurements. The residual between these two observations is significantly less than the residual between our measurements and the NVSS catalog. The additional discrepancies between our results and the NVSS database could result from a number of factors. For instance, a less accurate beam model or the flux variation of the NVSS source. Thus, the ~6.3% flux dispersion on average indicates an upper bound on the residual gain variations after the calibration process.

### 4.5.2. Flux Comparison with the Combined Map

Next, we make the flux density comparison using a map with a fine angular resolution created using the mapmaking process, as described in Haynes et al. (2018). We average the flux density of pixels within an aperture radius of $0.5\theta_{\text{FWHM}} \sim 1.5$

via (Fabello et al. 2011)

$$ S = \frac{\sum_i S(\theta_i)}{\sum_i K'(\theta_i)}, \quad \sigma_S^2 = \frac{\sum_i K'^2(\theta_i) \sum_i (S(\theta_i) - K'S)^2}{(\sum_i K'(\theta_i))^2 \sum_i K'^2(\theta_i)}, \quad (45) $$

where $\sigma_S$ is the measurements error, $\theta_i$ is the angular separation to the center pixel and $K'(\theta_i)$ represents the kernel function used in the map domain

$$ K'(\theta_i) = \exp \left[ -\frac{1}{2} \left( \frac{\theta_i}{\sigma_{K'}} \right)^2 \right], \quad (46) $$

where $\sigma_{K'} = \sigma_K + \theta_{\text{FWHM}}/(2\sqrt{2\ln 2})$, $\sigma_K$ is the kernel size of Equation (36) that applied during the mapmaking and $\theta_{\text{FWHM}}$ is the beam size. We use the same NVSS sources that were selected in Section 4.5.1 for comparison.

A direct comparison of bright pixels in the map with the NVSS source flux would show a large bias because the sources do not always transit across the beam center. However, no correction has been made in the mapmaking process because we cannot presume that we know the sources and their positions. With a sufficiently large number of scans, the sources would be completely sampled and the flux measurements taken after constructing the maps would be unbiased. However, due to the limited number of scans, and also RFI flagging, noise diode injection, and abandoning of data from bad beams, the sources are far from completely sampled. To recover the bias from the incomplete sampling of the sources, we simulate the TOD using the NVSS catalog. The simulated TOD has the same sky coordinates and mask as the real data, and is projected to the map domain using the same mapmaking procedure as used in the observations. We discover that the flux from most sources is pretty biased. The map-domain flux values are then corrected using the difference between the flux from the simulation and the NVSS

$$ S_{\text{FAST}} = S_{\text{map}} + (S_{\text{NVSS}} - S_{\text{sim}}), \quad (47) $$

where $S_{\text{map}}$ is the flux measured from the map, $S_{\text{NVSS}}$ is the flux read from the NVSS catalog, $S_{\text{sim}}$ is the flux measured from the map made by simulated TOD, and $S_{\text{FAST}}$ is the corrected flux. The comparison of the flux before and after correction is shown in Figure 25. The top panel shows the flux–flux comparison between the map-domain measurements and the NVSS catalog, and the bottom panel shows the histogram statistic of the relative flux residual defined in Equation (44). It is obvious that the measurements are significantly biased in the absence of flux correction. The flux correction makes a significant improvement, i.e., the measurement’s relative uncertainty is improved from 17.1% to 6.3% after the flux correction. The flux values of a total of 81 isolated sources are listed in the “$S_{\text{NVSS}}$” and “$S_{\text{FAST}}$” columns of Table 3. The correction fractions with respect to $S_{\text{NVSS}}$, i.e., $(S_{\text{NVSS}} - S_{\text{sim}})/S_{\text{NVSS}}$, are listed in the “Correction” column.

We also investigate how the aperture size affects the flux measurements. We vary the aperture radius size between 0.5 arcmin, 1.0 arcmin, and 1.5 arcmin, and show the comparison results in Figure 26. All of the measurements are corrected using the simulation with the corresponding aperture size. With varying aperture radius sizes, the flux values only slightly varied without a clear systematic trend.

The map-domain flux measurements result in about 6.3% relative uncertainty, which is consistent with the uncertainty of
<figure><img src="image_27.png" /><figcaption>**Figure 23.** The flux–flux diagram compares the measurements using the TOD with those in the NVSS catalog.</figcaption></figure>

the TOD flux measurements. This shows that our mapmaking procedure is accurate enough for continuum flux measurements. For the flux check in this work, we only selected a small number of bright and isolated point sources. We leave the work of identifying weak and diffuse sources to future studies.

**5. Summary**

Neutral hydrogen (HI) is known to trace galaxies in the post-reionization era. A comprehensive wide-field of the extragalactic HI survey could provide valuable information for both cosmology and astrophysics research. In this work, we report the TOD analysis pipeline designed for drift-scan observation using FAST.

The data analyzed in this work were collected over a few nights spanning 2019, 2020, and 2021. During the 4 hr drift scan of each night, the FAST telescope was pointed at a fixed altitude angle and the observation covered R.A. that ranged from 9 to 13 hr, which overlaps with the NGP area of SDSS. The FAST L-band 19-feed

<figure><img src="image_28.png" /><figcaption>**Figure 24.** The histogram of relative flux residual between our measurements using the TOD and the NVSS catalog, i.e., Equation (44).</figcaption></figure>

receiver is used in our observations. The feed array is rotated by 23°4 to optimize the coverage. The pointing decl. was shifted by 10'835 between different days to enlarge the survey area.

The noise diode signal, as the relative calibrator, is injected for 1 s in every 8 s. The noise diode signal is used to calibrate the bandpass gain and temporal drift of the gain. Our analysis indicates that the observation data in 2019 have significant bandpass shape variation during the 4 hr drift-scan observation. The bandpass shape of the data in 2021 is much more stable. The major data analysis focuses on the 7 nights of observations in 2021.

We applied the *SumThreshold* and *SIR* radio frequency interference (RFI) flagging program to the bandpass-calibrated data. To enhance the RFI signal and protect the potentially existing HI emission lines, the RFI flagging is applied to the feed averaged
<figure><img src="image_29.png" /><figcaption>**Figure 25.** Map-domain flux measurements comparing. The flux measurements before and after correction with simulation are shown in green and red, respectively. Top panel: The x-axis gives the flux values from the NVSS catalog and the y-axis gives the flux measurements from the map. Bottom panel: The histogram statistic of the relative flux residual, i.e., Equation (44). All of the measurements use an aperture radius size of 1.5 arcmin.</figcaption></figure>

data. Due to the contamination of the GNSSs, the data between the frequency range of .1135–1310 MHz are mostly flagged. In addition, about 40% and 14% of the data are flagged in the frequency range of 1050–1135 MHz (the low-frequency band) and 1310–1450 MHz (the high-frequency band), respectively.

We develop a temporal drift calibration strategy that estimates the gain variation across the drift-scan observation by applying a Wiener filter on the gain variation measurements. The Wiener filter is designed according to the 1/f noise temporal power spectrum model, which is constrained using the observation data. With our calibration strategy, a temporal oscillation of the gain is observed in the 4 hr drift scan and this oscillation can be well calibrated.

The absolute flux calibration is done by calibrating the noise diode spectrum against the celestial source 3C286. The calibration observation is made also in drift-scan mode. The noise diode spectrum shape is stable during the 7 nights observations in 2021. In addition, using the drift-scan observation of bright source 3C286, we check the beam profile for each of the feeds. The beam profiles of all of the 19 feeds are significantly asymmetric.

Due to the systematic background noise level variation during the observation time, a significant temporal baseline variation is observed with the gain-calibrated data. In particular, some sharp variations are shown in 1 night of the observations. These baseline variations are subtracted by fitting with a baseline template, which is constructed using the baseline averaged across different feeds. After baseline subtraction, the calibrated data are zero-centered and transferred to the standard mapmaking procedure.

<figure><img src="image_30.png" /><figcaption>**Figure 26.** Map-domain flux measurements with different aperture radius sizes. Top panel: The x-axis is the flux values from the NVSS catalog and the y-axis is the flux measurements from the map. Bottom panel: The histogram statistic of the relative flux residual, i.e., Equation (44).</figcaption></figure>

We check the standing-wave ripples across the frequency axis by estimating the delay spectrum using well-calibrated data. The standing-wave signature is more prominent in the bandpass measurements than in the sky observation, probably because this standing wave is induced by the noise diode. We use the low-pass filter with the size of $τ_c = 0.7$ µs to suppress both the noise and standing-wave signature during the bandpass determination. The standing-wave peak is negligible during our observations.

We check the measurement noise level using the TOD following the method introduced in Wang et al. (2021). The noise level of the low- and high-frequency bands are 116.39 mK and 99.47 mK, respectively. The noise level also varies between different feeds. The noise level of the feed in the outer circle of the feed array increases by more than 10% than the central feed. The noise level reduces slowly after integrating the measurements via mapmaking, due to weak RFI contamination, residual sky emission, or correlated noise.

We also study the systematic uncertainties by comparing the continuum flux measurements with the NVSS catalog. By applying the source-finding algorithm with the threshold of 30 µJy beam⁻¹, i.e., the flux limit due to the map rms, more than three thousand continuum sources are detected within our survey field. However, most of them are confused sources due to the angular resolution limit of the FAST beam. The flux-weighted differential number counts for the detected sources are consistent with the NVSS catalog down to ~7 mJy beam⁻¹, which is about five times the confusion limit. Finally, we chose 81 isolated NVSS sources with flux over 14 mJy beam⁻¹ in our survey field and find that the calibrated data shows about 4.2%, 5.3%, and 6.9%
measurement uncertainties for the central feed, inner-circle feeds, and outer-circle feeds, respectively. This uncertainty varies between the measurements of the central feed, inner-circle feeds, and outer-circle feeds. Finally, there is about 6.3% uncertainty on average, which is consistent with the map-domain flux measurements.

## Acknowledgments

This work made use of the data from the Five-hundred-meter Aperture Spherical radio Telescope (FAST). FAST is a Chinese national mega-science facility, operated by National Astronomical Observatories, Chinese Academy of Sciences. We acknowledge the support of the National SKA Program of China (Nos. 2022SKA0110100, 2022SKA0110200, 2022SKA0110203), the National Natural Science Foundation of

China (Nos. 11975072, 11835009, 12203061), the CAS Interdisciplinary Innovation Team (JCTD-2019-05), and the science research grants from the China Manned Space Project with No. CMS-CSST-2021-B01. L.W. is a UK Research and Innovation Future Leaders Fellow [grant MR/V026437/1].

## Data Availability

The data underlying this article will be shared on reasonable request to the corresponding author.

## Appendix

### Continuum Source Catalog

The flux density of the 81 isolated continuum sources over 14 mJy beam⁻¹ are list in Table 3.

<table><caption>**Table 3**
The Flux Density of the 81 Isolated Continuum Sources Over 14 mJy beam⁻¹ in the Survey Field</caption><thead><tr><th>RA (J2000)</th><th>Dec (J2000)</th><th>S<sub>NVSS</sub><br>[mJy]</th><th>S<sub>FAST</sub><br>[mJy]</th><th>Correction<br>[%]</th><th>RA (J2000)</th><th>Dec (J2000)</th><th>S<sub>NVSS</sub><br>[mJy]</th><th>S<sub>FAST</sub><br>[mJy]</th><th>Correction<br>[%]</th></tr></thead><tbody><tr><td>09h03m11s67</td><td>+26°09'40" 28</td><td>195.87</td><td>184.04 ± 3.88</td><td>13.35</td><td>09h42m36s39</td><td>+26°10'15" 27</td><td>35.09</td><td>27.41 ± 0.90</td><td>13.08</td></tr><tr><td>10h18m07s13</td><td>+26°09'42" 51</td><td>180.38</td><td>177.66 ± 3.80</td><td>13.99</td><td>09h26m00s53</td><td>+25°59'53" 30</td><td>113.56</td><td>115.60 ± 2.97</td><td>9.24</td></tr><tr><td>10h57m23s15</td><td>+26°01'30" 18</td><td>186.74</td><td>188.67 ± 4.89</td><td>4.92</td><td>11h34m02s19</td><td>+26°32'24" 54</td><td>15.62</td><td>13.75 ± 0.39</td><td>13.94</td></tr><tr><td>12h38m23s54</td><td>+26°33'46" 43</td><td>104.30</td><td>97.35 ± 3.34</td><td>-0.27</td><td>12h39m52s76</td><td>+26°32'46" 57</td><td>49.14</td><td>45.11 ± 1.52</td><td>10.42</td></tr><tr><td>12h48m41s35</td><td>+26°33'56" 48</td><td>74.42</td><td>70.56 ± 2.54</td><td>1.94</td><td>10h13m36s86</td><td>+26°52'06" 67</td><td>157.57</td><td>139.77 ± 3.66</td><td>9.56</td></tr><tr><td>10h45m57s36</td><td>+26°52'39" 93</td><td>56.87</td><td>46.13 ± 1.41</td><td>12.77</td><td>10h46m40s38</td><td>+26°54'15" 73</td><td>163.09</td><td>155.16 ± 3.46</td><td>13.96</td></tr><tr><td>09h21m58s57</td><td>+26°21'22" 60</td><td>206.52</td><td>211.16 ± 5.83</td><td>1.16</td><td>10h16m00s82</td><td>+26°22'00" 12</td><td>27.87</td><td>23.13 ± 0.81</td><td>2.94</td></tr><tr><td>10h48m28s59</td><td>+26°22'34" 89</td><td>237.69</td><td>244.67 ± 6.13</td><td>0.52</td><td>10h52m52s26</td><td>+26°22'16" 24</td><td>105.61</td><td>105.37 ± 2.62</td><td>2.66</td></tr><tr><td>09h27m45s09</td><td>+26°13'00" 51</td><td>150.04</td><td>146.74 ± 3.36</td><td>13.82</td><td>10h40m37s59</td><td>+26°13'24" 81</td><td>29.96</td><td>25.01 ± 1.14</td><td>13.45</td></tr><tr><td>10h54m56s64</td><td>+26°12'47" 26</td><td>338.51</td><td>337.11 ± 7.16</td><td>12.35</td><td>12h11m29s70</td><td>+26°14'50" 20</td><td>67.80</td><td>64.70 ± 1.25</td><td>14.18</td></tr><tr><td>09h49m06s44</td><td>+26°01'57" 71</td><td>20.78</td><td>11.63 ± 0.87</td><td>15.45</td><td>12h57m12s52</td><td>+26°01'40" 62</td><td>62.57</td><td>55.22 ± 1.23</td><td>1.79</td></tr><tr><td>09h18m56s17</td><td>+26°33'30" 99</td><td>20.90</td><td>18.22 ± 0.82</td><td>22.58</td><td>10h39m30s97</td><td>+26°33'03" 88</td><td>143.48</td><td>133.02 ± 3.53</td><td>1.31</td></tr><tr><td>10h59m07s98</td><td>+26°34'07" 50</td><td>56.83</td><td>55.56 ± 1.47</td><td>6.01</td><td>12h36m30s59</td><td>+26°35'19" 78</td><td>564.07</td><td>599.41 ± 20.89</td><td>29.70</td></tr><tr><td>11h22m04s96</td><td>+26°13'36" 73</td><td>303.49</td><td>300.70 ± 5.84</td><td>12.37</td><td>10h22m09s14</td><td>+26°23'13" 70</td><td>161.25</td><td>161.32 ± 4.29</td><td>7.96</td></tr><tr><td>10h55m47s35</td><td>+26°23'37" 82</td><td>152.02</td><td>150.63 ± 3.38</td><td>6.00</td><td>12h05m32s86</td><td>+26°25'40" 65</td><td>278.72</td><td>268.41 ± 4.18</td><td>27.54</td></tr><tr><td>11h25m00s85</td><td>+26°06'09" 32</td><td>56.43</td><td>46.94 ± 1.60</td><td>5.32</td><td>12h19m56s37</td><td>+26°07'04" 07</td><td>109.12</td><td>62.21 ± 2.72</td><td>10.48</td></tr><tr><td>12h21m02s00</td><td>+26°08'45" 74</td><td>199.08</td><td>123.17 ± 2.70</td><td>9.95</td><td>09h02m54s70</td><td>+25°55'44" 83</td><td>19.08</td><td>21.38 ± 1.02</td><td>19.49</td></tr><tr><td>09h16m47s38</td><td>+25°54'16" 05</td><td>73.79</td><td>68.79 ± 2.34</td><td>4.71</td><td>09h28m53s21</td><td>+25°54'37" 90</td><td>69.09</td><td>68.43 ± 2.45</td><td>9.84</td></tr><tr><td>09h33m15s67</td><td>+25°55'00" 69</td><td>61.57</td><td>62.32 ± 2.06</td><td>12.20</td><td>09h44m42s33</td><td>+25°54'43" 55</td><td>774.39</td><td>771.98 ± 19.53</td><td>7.81</td></tr><tr><td>09h45m53s28</td><td>+25°56'08" 37</td><td>54.64</td><td>55.52 ± 1.17</td><td>22.19</td><td>11h08m38s99</td><td>+25°56'13" 45</td><td>88.49</td><td>98.64 ± 2.68</td><td>5.22</td></tr><tr><td>09h40m13s47</td><td>+26°26'56" 68</td><td>124.22</td><td>107.93 ± 2.82</td><td>5.26</td><td>12h53m02s00</td><td>+26°27'45" 97</td><td>203.70</td><td>166.13 ± 7.02</td><td>2.76</td></tr><tr><td>09h18m37s49</td><td>+26°50'41" 10</td><td>105.95</td><td>103.85 ± 3.05</td><td>8.58</td><td>10h42m43s59</td><td>+26°50'02" 25</td><td>148.24</td><td>151.05 ± 3.48</td><td>6.64</td></tr><tr><td>12h04m12s95</td><td>+26°50'58" 95</td><td>417.72</td><td>420.74 ± 9.65</td><td>12.98</td><td>12h19m46s58</td><td>+26°17'34" 40</td><td>29.46</td><td>17.75 ± 0.92</td><td>17.36</td></tr><tr><td>12h20m27s95</td><td>+26°19'03" 54</td><td>32.41</td><td>17.92 ± 0.28</td><td>12.38</td><td>09h13m27s12</td><td>+26°02'23" 13</td><td>67.04</td><td>61.94 ± 1.10</td><td>30.81</td></tr><tr><td>09h17m15s57</td><td>+26°03'40" 53</td><td>58.42</td><td>49.17 ± 1.41</td><td>7.95</td><td>09h40m14s71</td><td>+26°03'29" 91</td><td>462.22</td><td>457.43 ± 13.37</td><td>21.65</td></tr><tr><td>09h41m28s66</td><td>+26°04'02" 56</td><td>79.37</td><td>77.97 ± 2.18</td><td>8.46</td><td>10h05m15s06</td><td>+26°03'35" 99</td><td>180.10</td><td>165.31 ± 3.55</td><td>23.18</td></tr><tr><td>10h24m19s44</td><td>+26°05'10" 71</td><td>20.86</td><td>18.40 ± 0.50</td><td>5.51</td><td>12h22m01s84</td><td>+26°04'15" 09</td><td>60.17</td><td>32.34 ± 0.98</td><td>28.94</td></tr><tr><td>09h30m27s63</td><td>+25°52'54" 26</td><td>78.58</td><td>83.92 ± 2.25</td><td>16.53</td><td>09h17m17s69</td><td>+26°45'51" 73</td><td>22.46</td><td>21.96 ± 0.73</td><td>14.87</td></tr><tr><td>10h35m16s54</td><td>+26°15'16" 45</td><td>388.16</td><td>367.07 ± 7.47</td><td>9.77</td><td>12h02m55s33</td><td>+26°15'18" 79</td><td>51.15</td><td>39.83 ± 1.00</td><td>13.88</td></tr><tr><td>11h34m34s12</td><td>+26°52'50" 95</td><td>108.72</td><td>99.49 ± 2.50</td><td>14.70</td><td>09h05m05s60</td><td>+26°11'23" 24</td><td>207.10</td><td>191.89 ± 4.47</td><td>10.57</td></tr><tr><td>10h11m06s75</td><td>+26°38'21" 01</td><td>33.34</td><td>27.02 ± 1.47</td><td>9.53</td><td>10h19m26s26</td><td>+26°38'58" 41</td><td>117.41</td><td>119.92 ± 3.15</td><td>10.87</td></tr><tr><td>10h25m08s32</td><td>+26°36'39" 24</td><td>14.53</td><td>8.41 ± 1.04</td><td>17.35</td><td>10h31m56s02</td><td>+26°38'20" 94</td><td>96.73</td><td>95.41 ± 3.55</td><td>9.19</td></tr><tr><td>10h57m24s93</td><td>+26°38'00" 02</td><td>435.24</td><td>420.13 ± 14.35</td><td>4.35</td><td>11h25m20s26</td><td>+26°38'02" 07</td><td>49.37</td><td>40.30 ± 1.67</td><td>5.63</td></tr><tr><td>09h23m43s10</td><td>+26°34'49" 94</td><td>41.15</td><td>38.81 ± 1.30</td><td>34.81</td><td>12h42m48s21</td><td>+25°58'40" 76</td><td>308.33</td><td>290.77 ± 6.32</td><td>17.88</td></tr><tr><td>11h06m08s35</td><td>+25°52'13" 04</td><td>38.26</td><td>39.91 ± 1.97</td><td>10.67</td><td>11h27m05s77</td><td>+25°50'55" 17</td><td>117.51</td><td>127.43 ± 3.32</td><td>2.96</td></tr><tr><td>12h00m05s07</td><td>+26°44'18" 34</td><td>34.49</td><td>40.04 ± 1.21</td><td>6.54</td><td>10h33m40s04</td><td>+26°42'57" 77</td><td>16.66</td><td>16.56 ± 0.57</td><td>13.50</td></tr><tr><td>10h54m40s25</td><td>+26°41'25" 69</td><td>63.90</td><td>59.84 ± 1.59</td><td>7.20</td><td>10h14m14s85</td><td>+25°52'00" 01</td><td>37.40</td><td>30.04 ± 1.06</td><td>21.30</td></tr><tr><td>09h54m39s79</td><td>+26°39'24" 55</td><td>312.67</td><td>302.08 ± 5.82</td><td>13.77</td><td>10h39m50s71</td><td>+27°01'48" 14</td><td>14.86</td><td>14.63 ± 0.73</td><td>5.01</td></tr><tr><td>12h24m21s36</td><td>+27°01'42" 20</td><td>97.89</td><td>95.13 ± 2.15</td><td>7.08</td><td>09h31m32s80</td><td>+26°40'48" 18</td><td>70.83</td><td>68.01 ± 2.04</td><td>12.48</td></tr><tr><td>11h02m53s30</td><td>+27°03'39" 38</td><td>35.19</td><td>28.59 ± 1.19</td><td>18.43</td><td>12h50m21s97</td><td>+27°04'32" 95</td><td>49.53</td><td>44.96 ± 1.43</td><td>18.40</td></tr><tr><td>09h15m01s99</td><td>+27°01'45" 69</td><td>26.15</td><td>30.06 ± 0.49</td><td>10.34</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

Note. The sources' sky coordinates at the J2000 epoch are presented in the "RA" and "Dec" columns; The flux density measurements read from the combined radio objects catalog (Kimball & Ivezić, 2008) are presented in the "S_NVSS" column; the flux density measurements from the map and corrected with the simulation are presented in the "S_FAST" column; and the "Correction" column presents the flux correction fractions with respect to S_NVSS, i.e., (S_NVSS - S Sim)/S_NVSS, where S Sim is the flux measured from the map made by the simulated TOD. The flux measurements use 1.5 arcmin aperture size and are corrected with the simulation.
ORCID iDs

Yichao Li https://orcid.org/0000-0003-1962-2013

Yougang Wang https://orcid.org/0000-0003-0631-568X

Furen Deng https://orcid.org/0000-0001-8075-0909

Wenxiu Yang https://orcid.org/0009-0006-2521-025X

Wenkai Hu https://orcid.org/0000-0002-3108-5591

Diyang Liu https://orcid.org/0009-0000-6895-9136

Xinyang Zhao https://orcid.org/0009-0008-2564-9398

Shifan Zuo https://orcid.org/0000-0003-3858-6361

Shuanghao Shu https://orcid.org/0009-0004-8919-7088

Jixia Li https://orcid.org/0000-0001-9652-1377

Peter Timbie https://orcid.org/0000-0003-0325-1633

Albert Stebbins https://orcid.org/0000-0002-3807-7252

Laura Wolz https://orcid.org/0000-0003-3334-3037

Fengquan Wu https://orcid.org/0000-0002-6174-8640

Xin Zhang https://orcid.org/0000-0002-6029-1933

Xuelei Chen https://orcid.org/0000-0001-6475-8863

References

Alam, S., Aubert, M., Avila, S., et al. 2021, PhRvD, 103, 083533

Amiri, M., Bandura, K., Chen, T., et al. 2023, ApJ, 947, 16

Anderson, C. J., Luciw, N. J., Li, Y. C., et al. 2018, MNRAS, 476, 3382

Anderson, L., Aubourg, É., Bailey, S., et al. 2014, MNRAS, 441, 24

Ansari, R., Campagne, J. E., Colom, P., et al. 2012, A&A, 540, A129

Bagla, J. S., Khandai, N., & Datta, K. K. 2010, MNRAS, 407, 567

Bandura, K., Addison, G. E., Amiri, M., et al. 2014, Proc. SPIE, 9145, 914522

Barnes, D. G., Staveley-Smith, L., de Blok, W. J. G., et al. 2001, MNRAS, 322, 486

Battye, R. A., Browne, I. W. A., Dickinson, C., et al. 2013, MNRAS, 434, 1239

Battye, R. A., Davies, R. D., & Weller, J. 2004, MNRAS, 355, 1339

Becker, R. H., White, R. L., & Edwards, A. L. 1991, ApJS, 75, 1

Bull, P., Ferreira, P. G., Patel, P., & Santos, M. G. 2015, ApJ, 803, 21

Chang, T.-C., Pen, U.-L., Bandura, K., & Peterson, J. B. 2010, Natur, 466, 463

Chang, T.-C., Pen, U.-L., Peterson, J. B., & McDonald, P. 2008, PhRvL, 100, 091303

Chen, X. 2012, IJMPS, 12, 256

Chen, Z., Chapman, E., Wolz, L., & Mazumder, A. 2023, MNRAS, 524, 3724

Cohen, A. S., Lane, W. M., Cotton, W. D., et al. 2007, AJ, 134, 1245

Cole, S., Percival, W. J., Peacock, J. A., et al. 2005, MNRAS, 362, 505

Colla, G., Fanti, C., Fanti, R., et al. 1972, A&AS, 7, 1

Condon, J. J. 1974, ApJ, 188, 279

Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693

Cunnington, S., Li, Y., Santos, M. G., et al. 2023, MNRAS, 518, 6262

Douglas, J. N., Bash, F. N., Bozyan, F. A., Torrence, G. W., & Wolfe, C. 1996, AJ, 111, 1945

Eisenstein, D. J., Zehavi, I., Hogg, D. W., et al. 2005, ApJ, 633, 560

Fabello, S., Catinella, B., Giovanelli, R., et al. 2011, MNRAS, 411, 993

Gao, L.-Y., Li, Y., Ni, S., & Zhang, X. 2022, arXiv:2212.08773

Giovanelli, R., Haynes, M. P., Kent, B. R., et al. 2005, AJ, 130, 2598

Giovanelli, R., Haynes, M. P., Kent, B. R., et al. 2007, AJ, 133, 2569

Gregory, P. C., Scott, W. K., Douglas, K., & Condon, J. J. 1996, ApJS, 103, 427

Harper, S. E., & Dickinson, C. 2018, MNRAS, 479, 2024

Harper, S. E., Dickinson, C., Battye, R. A., et al. 2018, MNRAS, 478, 2416

Haynes, M. P., Giovanelli, R., Kent, B. R., et al. 2018, ApJ, 861, 49

Hinton, S. R., Kazin, E., Davis, T. M., et al. 2017, MNRAS, 464, 4807

Hu, W., Li, Y., Wang, Y., et al. 2021, MNRAS, 508, 2897

Hu, W., Wang, X., Wu, F., et al. 2020, MNRAS, 493, 5854

Jarvis, M. J., Bhatnagar, S., Bruggen, M., et al. 2014, arXiv:1401.4018

Jiang, P., Tang, N.-Y., Hou, L.-G., et al. 2020, RAA, 20, 064

Jin, S.-J., Wang, L.-F., Wu, P.-J., Zhang, J.-F., & Zhang, X. 2021, PhRvD, 104, 103507

Kimball, A. E., & Ivezić, Ž. 2008, AJ, 136, 684

Lang, R. H., Boyce, P. J., Kilborn, V. A., et al. 2003, MNRAS, 342, 738

Li, D., & Pan, Z. 2016, RaSc, 51, 1060

Li, D., Wang, P., Qian, L., et al. 2018, IMMag, 19, 112

Li, J., Zuo, S., Wu, F., et al. 2020, SCPMA, 63, 129862

Li, Y., Santos, M. G., Grainge, K., Harper, S., & Wang, J. 2021, MNRAS, 501, 4344

Li, Y.-C., & Ma, Y.-Z. 2017, PhRvD, 96, 063525

Lidz, A., Furlanetto, S. R., Oh, S. P., et al. 2011, ApJ, 741, 70

Loeb, A., & Wyithe, J. S. B. 2008, PhRvL, 100, 161301

Mao, Y., Tegmark, M., McQuinn, M., Zaldarriaga, M., & Zahn, O. 2008, PhRvD, 78, 023529

Masui, K. W., Switzer, E. R., Banavar, N., et al. 2013, ApJL, 763, L20

Matshawule, S. D., Spinelli, M., Santos, M. G., & Ngobese, S. 2021, MNRAS, 506, 5075

McQuinn, M., Zahn, O., Zaldarriaga, M., Hernquist, L., & Furlanetto, S. R. 2006, ApJ, 653, 815

Meyer, M. J., Zwaan, M. A., Webster, R. L., et al. 2004, MNRAS, 350, 1195

Meyers, B. W., Hurley-Walker, N., Hancock, P. J., et al. 2017, PASA, 34, e013

Nan, R., Li, D., Jin, C., et al. 2011, IJMPD, 20, 989

Newburgh, L. B., Bandura, K., Bucher, M. A., et al. 2016, Proc. SPIE, 9906, 9906X5

Ni, S., Li, Y., Gao, L.-Y., & Zhang, X. 2022, ApJ, 934, 83

Offringa, A. R., de Bruyn, A. G., Biehl, M., et al. 2010, MNRAS, 405, 155

Offringa, A. R., van de Gronde, J. J., & Roerdink, J. B. T. M. 2012, A&A, 539, A95

Paul, S., Santos, M. G., Chen, Z., & Wolz, L. 2023, arXiv:2301.11943

Paul, S., Santos, M. G., Townsend, J., et al. 2021, MNRAS, 505, 2039

Perdereau, O., Ansari, R., Stebbins, A., et al. 2022, MNRAS, 517, 4637

Perley, R. A., & Butler, B. J. 2017, ApJS, 230, 7

Peterson, J. B., Aleksan, R., Ansari, R., et al. 2009, in astro2010: The Astronomy and Astrophysics Decadal Survey, 234

Pritchard, J. R., & Loeb, A. 2008, PhRvD, 78, 103511

Pritchard, J. R., & Loeb, A. 2012, RPPh, 75, 086901

Reid, B., Ho, S., Padmanabhan, N., et al. 2016, MNRAS, 455, 1553

Saintonge, A. 2007, AJ, 133, 2087

Santos, M., Bull, P., Alonso, D., et al. 2015, in Proc. Advancing Astrophysics with the Square Kilometre Array (AASKA14), 19

Santos, M. G., Bull, P., Camera, S., et al. 2016, in Proc. MeerKAT Science: On the Pathway to the SKA, 32

Seo, H.-J., Dodelson, S., Marriner, J., et al. 2010, ApJ, 721, 164

Spinelli, M., Carucci, I. P., Cunningham, S., et al. 2022, MNRAS, 509, 2048

Square Kilometre Array Cosmology Science Working Group, Bacon, D. J., Battye, R. A., et al. 2020, PASA, 37, e007

Sun, S., Li, J., Wu, F., et al. 2022, RAA, 22, 065020

Switzer, E. R., Masui, K. W., Bandura, K., et al. 2013, MNRAS, 434, L46

Tegmark, M. 1997, ApJL, 480, L87

Teunissen, P. J., & Montenbruck, O. 2017, in Handbook of Global Navigation Satellite Systems, ed. P. J. G. Teunissen & O. Montenbruck (Berlin: Springer)

Waldram, E. M., Yates, J. A., Riley, J. M., & Warner, P. J. 1996, MNRAS, 282, 779

Wang, J., Santos, M. G., Bull, P., et al. 2021, MNRAS, 505, 3698

Wolz, L., Blake, C., Abdalla, F. B., et al. 2017, MNRAS, 464, 4938

Wolz, L., Pourtsidou, A., Masui, K. W., et al. 2022, MNRAS, 510, 3495

Wu, F., Li, J., Zuo, S., et al. 2021, MNRAS, 506, 3455

Wu, P.-J., Li, Y., Zhang, J.-F., & Zhang, X. 2023a, SCPMA, 66, 270413

Wu, P.-J., Shao, Y., Jin, S.-J., & Zhang, X. 2023b, JCAP, 2023, 052

Wu, P.-J., & Zhang, X. 2022, JCAP, 01, 060

Wyithe, J. S. B., & Loeb, A. 2008, MNRAS, 383, 606

Wyithe, J. S. B., Loeb, A., & Geil, P. M. 2008, MNRAS, 383, 1195

Xu, Y., Wang, X., & Chen, X. 2015, ApJ, 798, 40

Zhang, M., Li, Y., Zhang, J.-F., & Zhang, X. 2023, MNRAS, 524, 2420

Zhang, M., Wang, B., Wu, P.-J., et al. 2021, ApJ, 918, 56

Zuo, S., Li, J., Li, Y., et al. 2021, A&C, 34, 100439

Zwaan, M. A., Meyer, M. J., Webster, R. L., et al. 2004, MNRAS, 350, 1210