\section{Introduction}
Over the past two decades, an increasing number of optical galaxy surveys, such as the 2DFGRS \citep{2001MNRAS.327.1297P}, SDSS\citep{1995wfsd.conf....3G,2004ApJ...606..702T,2009astro2010S.314S,2016AJ....151...44D}, WiggleZ\citep{2007ASPC..379...72G,2010MNRAS.401.1429D} are probing increasingly large volumes of the Universe, and provide large scale structure (LSS) data for cosmological studies. For example, assuming that the observed number density of galaxies traces the total density of the 
matter distribution, the baryon acoustic oscillation (BAO) features in the galaxy power spectrum are measured, and used as standard rulers to constrain cosmological models \citep{2005ApJ...633..560E,2005MNRAS.362..505C,2010MNRAS.401.2148P,2011MNRAS.416.3017B,2011MNRAS.418.1707B,2015MNRAS.449..835R,2017MNRAS.470.2617A,2018MNRAS.473.4773A}. While it is plausible that the galaxies formed from over-density perturbations and therefore trace the total density on large scales, it is vital to check this hypothesis and also understand the range of validity of it by observing the galaxies with different means. The 21cm line of the neutral hydrogen (HI) provides a good alternative way of observation in the radio wavelength. A number of HI galaxies surveys have also been carried out, e.g. the HIPASS survey \citep{2004MNRAS.350.1195M,2004MNRAS.350.1210Z}, and the ALFALFA survey \citep{2005AJ....130.2598G,2007AJ....133.2087S,2007AJ....133.2569G}, and JVLA deep survey \citet{2014arXiv1401.4018J}. However, limited by the  sensitivity of the telescopes, the redshift range of these surveys are much smaller than the current optical surveys.

A new generation of radio telescopes are being built or under development, including the Square Kilometer Array (SKA)  
in the southern hemisphere, and the FAST\citep{2011IJMPD..20..989N} in the 
northern hemisphere. These radio telescopes have much better sensitivities and observe HI galaxies at larger 
distances.  Here we consider FAST, which is about the complete its commissioning process and starts science runs.
 It has a very large aperture (300 meters during operation) and is to be equipped with multi-beam feed system and low-noise cryogenic
 receivers, ideal for conducting large surveys. 

The FAST has unprecedented large effective area and high sensitivity, nevertheless for a traditional 
galaxy survey \citep{2008MNRAS.383..150D} the redshift at which an individual galaxy could be detected is still very
limited, and its angular resolution would be insufficient to resolve the galaxy at high redshift. However,
to map the large-scale structure, in principle it is not necessary to resolve individual galaxies as 
traditional galaxy surveys do, instead the redshifted 21cm line intensity can 
be mapped with lower angular resolution, as is done in the Epoch of Reionization (EoR) 
experiments\footnote{One of us (XC) first realized that this mode of observation could be used with FAST to probe the large scale
structure and presented it at a meeting on the FAST science case held in May 2007 in Hangzhou, China.}.
In more general context, \citet{2008PhRvL.100i1303C} studied this mode of observation and named it the {\it intensity 
mapping} method, and  also proposed that  a cost-effective way to survey large scale structure is to 
develop a dedicated dense array of cylinder or small dish 
antennas \citep{2008PhRvL.100i1303C,2008arXiv0807.3614A,2010ApJ...721..164S,2012A&A...540A.129A}. Indeed
a number of such small-to-mid scale experiments are undergoing, such as those of 
Tianlai \citep{2012IJMPS..12..256C,2015ApJ...798...40X}, CHIME \citep{2014SPIE.9145E..22B} and 
HIRAX\citep{2016SPIE.9906E..5XN}, as well as the specially designed single dish experiment 
BINGO \citep{2012arXiv1209.1041B,2016arXiv161006826B}. 

For the FAST itself, several studies used the Fisher matrix formalism to make simple forecasts on the constraining power of 
cosmological parameters by HI galaxy survey \citep{2008MNRAS.383..150D} or intensity mapping surveys 
\citep{2015MNRAS.454.3240B,2017A&A...597A.136S,2019arXiv190803024Y}.

In this paper, we make a more detailed investigation by 
simulating the observed galaxies, and also compare the galaxy survey and intensity mapping. 
The layout of this paper is as follows. In Sec.II we describe our model of the telescope and its receiver feeds. the 
In section III, we present the modeling of the HI galaxies and their observation, as well as the simulated intensity map.  In section IV,  we make Fisher matrix
forecasts of the precision of power spectrum measurement using both HI galaxy surveys and intensity mapping surveys, and also make measurement using 
numerical simulation. The niche of HI galaxy survey and HI intensity mapping survey and the effect of foreground are discussed in Sec. V. Finally we summarize
the results in section VI.

\section{The FAST telescope}
In order to study how the FAST could survey the large-scale-structure, we conduct mock observations with 
simulated sky. We first generate a catalog of galaxies from simulation, then convert it into the simulated 
sky of HI intensity as would be observed by FAST. 

\subsection{The Instrument\label{sec:instr}}

The diameter of the FAST reflector is 500 m, the fully illuminated aperture at any time is $D=300$ m
since the telescope is designed to track 
objects. The beam size of the FAST is given by
\begin{eqnarray}
\theta &=& 1.22 \times \frac{21\ensuremath{\,{\rm cm}} (1+z)}{300 \ensuremath{\,{\rm m}}}= 2.94 (1+z) \arcmin	\end{eqnarray}
for observation of the 21cm line from redshift $z$.
 
During a drift scan, a single feed is fixed to be pointed to a particular declination in the due north or south direction, 
so that in a sidereal day, a ring of width $2.94(1+z) \arcmin$ centered at that declination is scanned. The pointing declination 
can be changed so as to cover the whole observable part of the sky. The FAST site is located at a latitude of 25$^{\circ}$48$^{\prime}$ North, 
and the maximum zenith angle is $40^\circ$,  allowing the observation of $\approx$ 50$\%$ of the full sky or about 20,000 $\deg^2$. 
  
 The FAST is equipped with a number of different feed and receiver systems. For HI survey, the most relevant 
are the L-band 19-beam feed/receiver system and 
the wide band receiver system. Additionally, there are also several low frequency receivers which cover down to 70 MHz,
which can be used for Epoch of Reionization (EoR) observations. Here we shall consider mainly the first two, which are relevant for 
low-or-mid redshift observations of large scale structure.  In addition, below we shall also consider a possible future UHF phased 
array feed (PAF) system. We summarize the information of these receiver systems in Table~\ref{receiver_table}. 

{\bf L-band 19 beam receiver system}. 
It covers the frequency range of 1.05-1.45 GHz, and the beams are arranged in two 
concentric hexagonal rings around the central beam. 
The minimum spacing between beam centers is 5.73 arcmin and is approximately constant, though for each beam the width scales roughly as $\theta \propto (1+z)$.
In this paper we assume the feed array are tilted an angle of 23.4 degree
with respect to the compass points to increase the area covered for each scan, 
as was proposed for the Commensal Radio Astronomy FasT survey (CRAFTS)\citep{2018IMMag..19..112L}, though we note that this is not the only choice available.
The whole 19 beams span 22.8 arcmins across the north-south direction at 1.42GHz (calculated for the centre of the beam).  The sky is covered by shifting the whole array in declination by 21.9 arcmin for the next scan.
A drift scan of $\pm 40^\circ$ from the center declination would require about 220 strips 
(i.e. 220 days)  to cover the region once. 



\begin{table}
	\centering
	\caption{FAST survey receiver parameters. The $t_{\rm sur}$ refers to the time needed to finish a full drift scan of $\pm 40^\circ$ from 
	the center declination of FAST. The $T_{\rm rec}$ is the receiver noise.}
	\label{receiver_table}
	\begin{tabular}{ccccc} 
		\hline
		receiver & band(GHz) &  Beams &$T_{\rm rec}$(K)& $t_{\rm sur}$(days)\\
		\hline
		L-band & 1.05-1.45 & 19 & 20 & 220\\
        Wide-band & 0.27-1.62 & 1 & 60 & 1211\\
        UHF PAF (future) & 0.5-1.0 & 81 & 30 & 135 \\
		\hline
	\end{tabular}
\end{table}

{\bf Wide-band receiver system.} 
For higher redshift (z>0.35), at present the survey can be done with a single feed wide band receiver, which covers a frequency 
from 0.27GHz to 1.62 GHz.The receiver noise for Wide-band receiver system is $\approx$ 60K, to have the same noise scale, the survey for redshift larger than 0.35 need twice more time than the survey for redshift smaller than 0.35.  The strip width in this case is $2.9(1+z)\arcmin$, so to 
cover the $\pm 40 \deg$ sky, it would require 1211 days with single feed strip to cover the same sky region at $z=0.35$, which is much less practical due to the 
long observation time required. 

{\bf The PAF receiver system for UHF band.}
 In the future, it is worthwhile to consider equipping the FAST telescope 
with a multi-beam receiver at the lower frequency band for a survey of higher redshifts. 
A phased array feeds (PAF) with cryogenic receiver system would allow rapid survey of large areas of sky, and such 
development has also been pursued for FAST  \citep{wu2016}. 
Here as an illustrative example, we consider a low frequency PAF system with 81 effective beams, 500MHz bandwidth centered at 0.75 GHz,  
(i.e. 0.5-1.0 GHz), a system noise about 30K,  and an aperture efficiency around 70$\%$. These beams could positioned in square array, then a pixel in the sky 
will be scanned 9 times if drift along a side of the square. In this way, a full drift scan of $\pm$40$^{\circ}$ could be finished in 135 days, 
with integration time of 291s on each pixel. 


\subsection{Integration Time and Noise}
The sky drifts across with a speed of $\omega_e \cos \delta$ in a drift scan survey,  where 
$\omega_e \approx 0.25 \arcmin/s$ is the angular velocity of the rotation of the
Earth, and $\delta$ is the declination of the pointing.  The time for drifting across a pixel is given by 
\begin{eqnarray}
  t_{\rm pix}&=& 2.9(1+z) \arcmin /(\omega_e \cos\delta).
\end{eqnarray}
One circle is completed in a sidereal day, though in practice the night time data is usually of 
much smaller noise than the day time data.  
At $z=0$, and $\cos\delta \approx 1/2$ (near the zenith of FAST site), we get 24s per beam. And because of the overlap of 19 beams in one horizontal scanning(see Fig. \ref{noise_strip}), most pixels in a 19-beams strip will be scanned twice, resulting in 48s per beam. Within the observable part of sky, the circles with 
higher declination (northern part of sky) has smaller area, while the integration time per pixel is larger. The expected thermal noise for a dual polarization single beam is
\begin{eqnarray}
\sigma_{\rm noise} = \sqrt{2}\frac{k_B T_{\rm sys}}{A_{\rm eff}}\frac{1}{\sqrt{\Delta\nu t}}
\label{thermal_noise}
\end{eqnarray}
where $t$ is the total integration time, $\Delta\nu$ is the frequency bandwidth for a channel, and $k_B$ the Boltzmann constant. The aperture efficiency is about 70$\%$, giving effective aperture $A_{\rm eff} \approx$ 50000 $\ensuremath{\,{\rm m}}^2$. 
The system temperature is $T_{\rm sys}= T_{\rm rec} + T_{\rm sky}$, 
where $T_{\rm rec}$ is the receiver temperature and taken to be 20K, and away from the Galactic plane, the sky temperature is modeled as 
\begin{eqnarray}
T_{\rm sky} = 2.73 + 25.2 \times(0.408/\nu_{\rm GHz})^{2.75} K.
\label{sky_noise}
\end{eqnarray}

 If we assume a velocity line width of $\approx$ 5 $\kms$ for the 
spectral line observation in HI galaxy, and 48s integration time per beam, the instantaneous sensitivity of each beam of the FAST system will be 0.86mJy. Below we shall consider surveys with an average of 48s, 96, 192s and 384s integration time per beam, according to once, twice, three times and  four times repeat observation respectively.

In the case of 19 beam L-band feed,  approximately every pixel would be covered by several beams, effectively double the integration time. In a more careful treatment,  we may estimate the noise as follows. The time stream data is related to the signal by 
\begin{eqnarray}
\mathbf{d={A s} + n},
\end{eqnarray}
where the time-ordered data vector $\mathbf{d}$ has a dimension of $19 N_t$ 
where $N_t$ is the length of the time-ordered-data, the sky pixel vector $\mathbf{s}$ has a dimension $N_{\rm pix}$, 
and the pointing matrix $\mathbf{A}$ has a dimension of $19 N_t \times N_{\rm pix}$. 
The minimum variance estimator for the sky is 
\begin{eqnarray}
\mathbf{ \hat{s} = (A^t N^{-1} A)^{-1} A^{t} N^{-1} s}
\end{eqnarray}
where $\mathbf{N}$ is the covariance matrix of the noise in the time-ordered data.
The sky map noise covariance matrix is then 
\begin{eqnarray}
    \mathbf{C_{N} = (A^{T}N^{-1}A)^{-1}}.
\label{noise_matrix}
\end{eqnarray}
Using this expression we can estimate the map noise.  

\begin{figure}
  \centering
     \includegraphics[width=0.45\textwidth]{noise_strip.eps}
     \includegraphics[width=0.45\textwidth]{noise_strip_repeat.eps}
  \caption{The top panel shows the strips of noise level produced by one scan of the 19-beam feed. The bottom panel shows the noise level from multiple scans with vertical intervals of 21.9 arcmin. The width of strip is selected to be 2.3 arcmin, which is the vertical intervals of two nearest dashed line in the top panel.}
  \label{noise_strip}
\end{figure}

In Fig. \ref{noise_strip}, we show the estimated noise of the sky map obtained by the 19-beam receiver in units of 
single beam receiver. For simplicity we assumed the beams are identical and have a gaussian beam within the beam width, 
though in reality there is much difference in the central and outer beams. Also, we assumed a constant system temperature, though actually the system temperature varies, as the sky temperature varies. 

As one might expect, in a scan along the horizontal direction of the 
19 beam receiver, pixels which are near the center of the receivers will be scanned by more than one beams, 
resulting in a lower noise than others. Such inhomogeneous noise distribution is undesirable, because in the large scale structure
measurement it may bias the observation and induce superfluous structures. The bias may be approximately corrected by
introducing selection functions, but as the real noise is varying and not accurately known, precision is hard to achieve. To reduce 
such effects, we need to have a relatively uniform distribution of noise in the survey regions. From Fig. \ref{noise_strip},
we see by partially overlap the scanning strips (with a vertical intervals of 21.9 arcmin), a large part of this inhomogeneity 
could be removed, making it a nearly uniform survey in the central part.
\section{Simulation}

\subsection{The Galaxy Model}
We used the catalog from the Semi-Analytic Suite of the SKA Simulated Skies($S^3$-SAX), in which 
the cosmic evolution of the galaxies is tracked by semi-analytic models \citep{2007MNRAS.375....2D} based on the 
Millennium N-body simulation \citep{2005Natur.435..629S}, and the amount of neutral atomic hydrogen (HI)
and molecular (H$_{2}$) hydrogen in galaxies are computed with the semi-analytical 
model \citep{2009ApJ...698.1467O,2009ApJ...702.1321O,2009ApJ...703.1890O}. 
An easy-to-use mock catalog (\citet{2014arXiv1406.0966O}) of galaxies with detailed physical 
properties (position in the sky, apparent redshift, stellar mass, HI mass, effective radius, etc.) is available. The catalog is for a cone 
with a field of 10-by-10 degrees and a redshift range of 0.0-1.2. 
It is complete down to an HI mass of $10^8 \ensuremath{\, {\rm M}_{\odot}}$. A deficiency of this model is that this HI cutoff mass is still relatively high, which could miss a significant amount of HI in dwarf galaxies. 
This limit is only a minor concern when dealing with isolated direct HI detections in blind 
surveys, because only a tiny fraction of the total survey volume is sensitive to HI masses < $10^8 M_{\sun}$. 
However, when dealing with global HI mass estimates, e.g. the intensity mapping experiment, the 
HI mass contained in unresolved galaxies is non-negligible. 
We compute $\Omega_{\rm HI}(z)$ from all the galaxies in the mock catalogue, and compare it with the 
observations \citep{2005MNRAS.359L..30Z,2007MNRAS.376.1357L,2010ApJ...723.1359M,2011ApJ...727...40F,2012ApJ...749...87B,
2013MNRAS.433.1398D,2006ApJ...636..610R,2017MNRAS.471.3428R,2013MNRAS.435.2693R,2016MNRAS.460.2675R,
2018MNRAS.473.1879R,2015MNRAS.452.3726H,2016ApJ...818L..28K,2016ApJ...818..113N,2018MNRAS.tmp..502J,2019MNRAS.489.1619H},  
the result is shown in  Fig. \ref{omega_fit}. There are still quite large scatters and discrepancy in the result, but already we can see
the $S^3$-SAX simulation may have under-estimated the amount of HI by a factor of between 1.2 to 2.0,
especially at higher redshifts. The more recent MUFASA cosmological hydrodynamical simulation \citep{2017MNRAS.467..115D}
is in better agreement 
with the observations. We have done most of our galaxy survey simulation with the $S^3$-SAX mock catalogue, for 
the computation of the intensity mapping, we scale all HI flux(equivalent to HI mass) with a 
z-dependent factor to compensate the lost HI mass.

\begin{figure}
  \centering
  \includegraphics[width=0.45\textwidth]{omega_fit.eps}
  \caption{The $\Omega_{\rm HI}(z)$ from the mock catalog (blue line, and red line for the smoothed result), and the 
  various observations in . The dashed green line represents results in 
  the MUFASA cosmological hydrodynamical simulation from: $\Omega_{\rm HI}= 10^{-3.45}\times(1+z)^{0.74}$.}
  \label{omega_fit}
\end{figure}


For each galaxy in the catalog, the HI distribution and 21cm emission is modeled following \citep{2016MNRAS.460.4366E}. 
We generate a mini data cube for each galaxy, then we re-grid the mini data cube to a full-size cube which contains 
all the galaxies . 

The HI mass density distribution of the galaxy is modeled as a thin axisymmetric exponential model:
\begin{eqnarray} 
\Sigma_{\rm HI}(r) = \frac{\tilde{\Sigma}_{\rm H} e^{-r/r_{\rm disk}}}{1+R^{c}_{\rm mol} e^{-1.6r/r_{\rm disk}}},
\label{sigma_HI}
\end{eqnarray}
where $r$ denotes the galactocentric radius, $r_{\rm disk}$ refers to the scale length, $\Sigma_{\rm HI}(r)$ is the surface 
density of the total hydrogen component, $\Sigma_{\rm H} = M_{\rm H} /(2\pi r^{2}_{\rm disk} )$ is a normalization factor 
and $R^{c}_{\rm mol}$ denotes the ${\rm H}_2$ /HI-mass ratio at the galaxy center. 
This is based on a list of empirically supported assumptions, (i) the cold gas of regular galaxies resides in a flat disk 
(see \cite{2008AJ....136.2782L} for local spiral galaxies, \cite{2002AJ....124..788Y} for local elliptical galaxies, 
\cite{2006ApJ...640..228T} for galaxies at higher redshifts); (ii) the surface density of the total hydrogen component (HI+H2) 
is well described by an axisymmetric exponential profile (\cite{2008AJ....136.2782L}); (iii) the local H2/HI-mass ratio scales 
as a power of the gas pressure of the ISM outside molecular clouds (\cite{2006ApJ...650..933B}). Recent 
observation(\cite{2012MNRAS.422.1835S,2017MNRAS.464.1029N}) for early type 
galaxies reveals that for most of the elliptical galaxies, the HI gas has a morphology that is similar in appearance to the 
discs of radio emission associated with SF in spiral galaxies. Most of the HI detection exhibits a large, settled HI disc or ring.
The orientation of the galaxy is randomly chosen. In reality, the galaxy may have some correlation of 
intrinsic alignment in their orientation, but such alignment is generally a second order effect, and does not significantly 
affect the analysis given below.
When generating a galaxy model with the parameters from the catalog, we convert the apparent HI half-mass radius, 
$R_{\rm HI}^{\rm half}$, along the major axis into an exponential disc scale length $r_{\rm disk}$. The galaxy is modeled 
out to a radius of 3.5 $r_{\rm disk}$. With the surface density deduced, we then convert it to the mass distribution. 
The circular velocity profile of the galaxy is modeled with the Polyex analytic function \citep{2002ApJ...571L.107G}: 
\begin{eqnarray}
V_{\rm PE}(r) = V_{0}(1-e^{-r/r_{\rm PE}})(1+\ \frac{\alpha r}{r_{\rm PE}})
\label{circular}
\end{eqnarray}
where $V_{0}$, $r_{\rm PE}$, and $\alpha$  determine the amplitude, exponential scale of the inner region, 
and the slope of the outer part of the rotation curve respectively. These parameters are derived from the luminosity of the 
galaxy given in the semi-analytical model,  using the empirical relations derived from nearly 2200 low redshift disk galaxies\citep{2006ApJ...640..751C}.
The semi-analytical model gives R-band luminosity, while the \citet{2006ApJ...640..751C} model used I-band luminosity, 
so we convert them by $M_{I} = M_{R} - 0.37$ \citep{2012MNRAS.426.3385D}.



The HI flux density is then given by
\begin{eqnarray}
\frac{M_{\rm HI}}{M_{\sun}}= 2.36 \times 10^5 (\frac{D_L}{\ensuremath{\,{\rm Mpc}}})^2 \frac{S_{i}}{\ensuremath{\,{\rm Jy}}} \frac{d v}{\kms} (1 + z)^{-2},
\label{HI_flux}
\end{eqnarray}
where $S_{i}$ is flux density in units of Jy in channel i of the mini-cube, dv is the velocity width of a channel in$\kms$, and $D_{L}$ is the luminosity distance of the target galaxy in $\ensuremath{\,{\rm Mpc}}$ units, and z is its evaluated redshift. Note that $dv$ here is defined in intrinsic velocity bin, 
if it is for the observed velocity bin, the $(1 + z)^{-2}$ will be replaced by $(1 + z)^{-1}$. 
For each galaxy, the  data cube has $100 \times 100$ angular pixels, and $5\kms$ in velocity channel width. Each voxel (volume pixel) have its HI flux density and velocity that computed, and we reposition it in a 3D data cube.


We use a light cone catalogue from the \citet{2014arXiv1406.0966O} simulation, which spans a field of $10\times10 \deg^2$ on the 
sky and a redshift range of 0.0-1.2. This volume contains 19,210,309 galaxies with a total HI mass 
of $2.065 \times 10^{16} M_{\sun}$. 
We then re-grid the mini data cubes into the full-size synthetic cube and place it in the corresponding angular position and frequency. The 
final full-size synthetic cube have a pixel width of 0.0133 deg 
and a fixed channel width of 
0.0237MHz, corresponding to 5$\kms$ at the redshift of z = 0.
To simulate the sky observed by FAST, we convolve each channel of the synthetic data cube with a redshift-dependent Gaussian 
point-spread-function (PSF), with beam width proportional to $1+z$. 




\subsection{HI Galaxy Detection}
To simulate the detection of HI galaxies, we first re-bin the full-size synthetic data in RA and DEC to 
an angular resolution of 0.08degree, corresponding to 2 times FWHM of the FAST beam. 
If we assume the noise for each beam is gaussian, the 
noise for each pixel is rescaled as 
\begin{eqnarray}
    \sigma_{\rm noise}^{\rm pixel} = \sigma_{\rm noise}^{\rm beam}\times\sqrt{\frac{A_{\rm pixel}}{A_{\rm beam}}},
\label{thermal_noise_pixel}
\end{eqnarray}
where $A_{\rm pixel}$ is the sky area of a pixel and $A_{\rm beam}$ is the sky area of the beam.

In the frequency axis, the synthetic data cube is re-binned to a resolution of 0.1MHz, 
corresponding to a velocity width of 20km/s at z = 0. Because the flux and noise scale will change with 
the bandwidth, if the full-size synthetic data is smoothed to a velocity width of $W_{s}$$\kms$, 
the signal-to-noise from one velocity bin is scaled as (\citet{2007AJ....133.2087S})
\begin{eqnarray}
    S/N = \frac{(F_{0}/N)/W_{s}}{\sigma_{0}}\times\bigg(\frac{W_{s}}{5\ensuremath{\,{\rm km}} \mathrm{s}^{-1}}\bigg)^{1/2},
\label{signal_to_noise}
\end{eqnarray}
where $F_{0}$  is the total velocity integrated HI flux of a galaxy, $N$ is the number of the velocity bins the galaxy spans 
and $\sigma_{0}$ is the thermal noise with a velocity width of 5$\kms$. 

With the above mock data, we may simulate galaxy detection as follows:
\begin{enumerate}
 \item {\it coarse resolution search.} Re-bin the noise filled data to a angular resolution of 0.08$\deg$(two times of FWHM of FAST) 
 and a frequency resolution of 0.473MHz (corresponding to velocity resolution of 100km/s at redshift 0), 
 setting a threshold of $3\sigma$, and detect voxels above the threshold.
 \item {\it fine resolution fit.} For galaxies detected in the coarse search, use a finer frequency resolution (0.0236MHz, corresponding to 
 velocity resolution of 5km/s at redshift 0) to fit its spectrum in the data cube with a parameterized profile function. 
 If a reasonable HI profile is obtained,  we integrate the HI profile, the candidate is selected as a galaxy if the total flux exceeds 5 $\sigma$.
 About  20$\%$ of the candidates found in the first step passed the second step, the other ones might be large noise.
 \end{enumerate} 


For the HI profile, we use the so-called "busy function" proposed by  \citet{2014MNRAS.438.1176W}, 
which has great flexibility in 
fitting a wide range of HI profiles from the Gaussian profiles of dwarf galaxies to the broad, asymmetric double-horn profiles of 
spiral galaxies. 



\subsection{Galaxy Distribution}
To show the capacity of surveys with different integration time, we have run our selection pipeline with three different $\sigma_{0}$: 0.86mJy, 0.61mJy, 0.43mJy and 0.31mJy, corresponding to 48s, 96s, 192s and 384s integration time per beam.

 \begin{figure}
  \centering
  \includegraphics[width=0.45\textwidth]{polar_distribution_10Mpch.eps}\\
  \includegraphics[width=0.45\textwidth]{polar_distribution_all.eps}
  \caption{The top panel is the 10 $\ensuremath{\,{\rm Mpc}}/h$ slab of the galaxy candidates from the FAST 2-years(or 192s-integration-time) galaxy survey, 
  the bottom panel is the distribution of all galaxies in the catalog.}
  \label{polar_distribution}
\end{figure}

We plot the 10 $\ensuremath{\,{\rm Mpc}}/h$ slab of the mock galaxies detected in a two year (or 192s-integration-time) survey in 
Fig.~\ref{polar_distribution}, all galaxies in the simulated are also plotted for comparison. The masses and redshift of these detected galaxies are plotted in Fig.~\ref{mass_dist_distribution_multitime}, and in Fig.~\ref{number_density_multitime} we show 
the number density of the detected galaxies.
We can see from these figures that as redshift increases, the galaxies thinned out, and 
the number density decrease drastically as redshift increases. At $z=0.2$, the number density already falls off by two or three orders of magnitude compared with the nearby galaxies, and only the most 
massive galaxies can be detected at the higher redshifts.  Indeed, 
existing HI galaxy surveys are all limited to $z<0.2$, and even FAST, 
the largest single dish telescope in the world,  could not detect many galaxies in this mode. 
In the case of integration time of 48 s, the galaxy number density would drop to about $10^{-4} (\ensuremath{\,{\rm Mpc}}/h)^{-3}$ at $z=0.15$, 
corresponding to roughly the galaxy number density required for BAO measurement.
For the longer integration time, the distribution of the detectable galaxies extends to higher redshift. 
But even for the 384 s survey, the number density drops to $10^{-4} (\ensuremath{\,{\rm Mpc}}/h)^{-3}$ at $z=0.4$.

\begin{figure}
  \centering
  \includegraphics[width=0.45\textwidth]{mass_z_distribution_volumelimited.eps}
  \caption{The distribution of the HI mass of the galaxy candidates from the 
  15$\times$15$\times$600 $h^{-3}\ensuremath{\,{\rm Mpc}}^3$ comoving volume.}
  \label{mass_dist_distribution_multitime}
\end{figure}


\begin{figure}
  \centering
  \includegraphics[width=0.4\textwidth]{number_density_multitime.eps}
  \caption{The comoving number density computed from the detected catalog.}
  \label{number_density_multitime}
\end{figure}



\begin{figure}
  \centering
  \includegraphics[width=0.4\textwidth]{confusion_rate_multitime.eps}
  \caption{The confusion rate of the detected galaxies from different surveys. The minimum HI mass is set as $10^{8.5}M_{\odot}$.}
  \label{confusion_rate_ratio_multitime}
\end{figure}

For galaxy survey, galaxies may overlap with each other along the line of sight. The ability to uniquely identify the 
individual galaxy is an important evaluation of the performance of a telescope. To quantify this, we introduce the 
confusion rate, which is defined as the fraction of galaxies fall within the same voxel.
The voxel have a length of $0.08 \deg$ a side, which is 2 times of the FWHM of the FAST beam at $z=0$ , and a 
bandwidth of 1MHz in frequency axis. We have set a limit of HI mass of $10^{8.5}M_{\odot}$ in the 
mock catalog to remove the effect induced by low mass galaxies. 
For the galaxy detection algorithm  we used,  the confusion rate is shown 
in Fig. \ref{confusion_rate_ratio_multitime}. As the redshift increases, the confusion rate rises rapidly.
The integration time has a mild influence on the confusion rate, as the smaller galaxies are detected for longer
integration time, which leads to an increase in the number of galaxies per beam and hence a higher confusion rate.
This steady raise of the confusion rate shows how the galaxy count transits smoothly and naturally to 21cm intensity 
mapping with the increase of redshift in the FAST observation. 

For the detected galaxies, due to the finite beam size and the noise, there will also be errors in the measured positions.
To quantify this effect we compare the measured positions in the mock observation with the original position in the 
catalog. The measured central position is taken to be the flux weighted average, 
\begin{eqnarray}
    \nu_{c} = \frac{\int s(\nu)\nu d\nu}{\int s(\nu)d\nu}, \qquad \theta_c = \frac{\int s_{i}\theta_{i} d\theta}{\int s_{i}d\theta},
    \label{freq_c}
\end{eqnarray}
The frequency integration range covers about $\pm$200 km/s, with resolution of 0.0236MHz.  
We sum the flux from the sub-data cube along the frequency axis, the angular beam average is computed 
on a grid with spacing of 0.0133$\deg$, and beam width 0.08$\deg$. 

The error of the positions measurement in comoving coordinates is shown in Fig. \ref{position_shift}. We see that the noise 
in the galaxy HI profile can induce a shift in its position. For the plotted redshift range ($z<0.3$), the shift in the direction perpendicular to the line of sight (los) is generally smaller than that along the los direction. About 95$\%$ galaxies have position shifts in the perpendicular direction lower than 0.1$\ensuremath{\,{\rm Mpc}}/h$ and about 90$\%$ galaxies have shift in parallel direction lower than 0.5$\ensuremath{\,{\rm Mpc}}/h$. We colored the galaxies with their velocity-integrated flux. It shows that the galaxies with low velocity-integrated flux tend to have large shift. There are more galaxies which have large shift in survey with higher integration time, because more galaxies with low velocity-integrated flux can be detected with longer integration time. 
We also show the standard deviation of the position shift with blue vertical lines at different redshifts in 
Fig. \ref{position_shift}. Higher redshift has larger shift because there are more galaxies with lower velocity-integrated flux.

\begin{figure*}
      \centering
      \includegraphics[width=0.7\textwidth]{delta_xy.eps}
      \includegraphics[width=0.7\textwidth]{delta_z.eps}
      \caption{The position error induced by thermal noise in the observation for different integration times. Top panels: the 
      error perpendicular to the los. Bottom panels: the error along the los.  The galaxies are 
      colored according to their integrated flux.}
      \label{position_shift}
\end{figure*}
\begin{figure*}
  \centering
  \includegraphics[width=0.8\textwidth]{flux_with_noise_192s.pdf}
  \caption{Noise-filled HI flux slice from different redshifts, all have the same frequency interval 1 MHz. The galaxy flux and noise both are computed with a bandwidth of 1 MHz and a integration time of 192s per beam.}
  \label{HI_intensity_flux_noise}
\end{figure*}
\begin{figure*}
  \centering
  \includegraphics[width=0.8\textwidth]{HI_intensity.pdf}
  \caption{Velocity integrated HI flux map from different redshifts with $\Delta z = 0.1$.}
  \label{HI_intensity}
\end{figure*}


\subsection{HI Intensity Map}
In Fig.~\ref{HI_intensity_flux_noise} we show slices of HI intensity maps at six different redshifts in the range $0<z<0.3$, the depth is 1 MHz for each slice. We also added a thermal noise corresponding to 1 MHz of channel bandwidth 
and 192 seconds of integration time. To obtain the noise-filled maps, we construct the flux cube and noise cube with a bandwidth of 1 MHz, respectively, then combine them together. As we have mentioned above, missing the galaxies with HI mass below  $10^8 M_{\sun}$ will have a significant 
influence in intensity mapping experiment. We correct this by scaling all HI masses in the mock catalog with a z-dependent
 factor (> 1) to match $\Omega_{\rm HI}$ inferred from MUFASA cosmological hydrodynamical simulation. 
As can be seen from this figure, even for the relatively long integration time of 192 seconds, 
the map is still largely dominated by noise,  showing the challenges one will face in HI surveys.

In Fig.~\ref{HI_intensity}  we show the projected HI intensity maps at $0<z<0.8$, each with a redshift interval of 0.1.
The HI distribution is shown more clearly in these maps without the noise. For slices nearby, one can clearly see individual galaxies. As the distance increases, the structures become more 
blurred, and also the intensity drops. Individual galaxies become increasingly
difficult to see, but the overall structure remains, which illustrates how the intensity mapping could be used to probe the large scale structure. Note that the equal spacing in redshift means
slightly smaller comoving distance spacing at higher redshift, but we have checked and found that the equal spacing 
in comoving distance generate maps pretty similar to these.


\section{Power Spectrum Measurement}
The power spectrum is the most widely used statistics for large scale structure. In this section we describe its measurement and 
error forecast for both the HI galaxy survey and the HI intensity mapping survey with FAST.

\subsection{HI Galaxy Power Spectrum}

In a galaxy redshift survey with negligible error on the position of galaxies, 
the measurement error of the power spectrum comes from sample variance as well as shot noise.
Over a k bin of width of $\Delta k$ is \citep{1994ApJ..426...23F,2008MNRAS.383..150D}
\begin{eqnarray}
\frac{\sigma_{P}}{P} = \sqrt{2\frac{(2\pi)^3}{V_{\rm eff}}\frac{1}{4\pi k^2 \Delta k}}\frac{P(k)+1/n}{P(k)},
\label{projectederror_galaxy}
\end{eqnarray}
where
$V_{\rm eff}(k) = \int{\bigg[\frac{n(\vec{r})P(k)}{n(\vec{r})P(k) + 1}\bigg]^2}d^3\vec{r}$ for detected galaxies $n(r)$.
The error on galaxy position may also induce slight errors, but for the scale of interest they are negligible.

\begin{figure}
  \centering
  \includegraphics[width=0.4\textwidth]{projected_error_galaxy_count.pdf}
  \caption{The projected error on power spectrum from a $20000\deg^2$ galaxy survey with integration time of 48 s, 96 s, 192 s and 384 s per beam. 
  At k$\approx$ 0.07 $h/\ensuremath{\,{\rm Mpc}}$ the S/N can reach 5.0 at $z\approx0.2, 0.25, 0.3$ and 0.35 respectively.} 
  \label{projected_error_galaxy_count}
\end{figure}

An optimal weighted estimator (known as the Feldman-Kaiser-Peacock or FKP estimator)
may be formed to minimize the measurement error of the power spectrum \citep{1994ApJ..426...23F}.
To make the measurement  in the irregular geometry of an actual survey, a mock sample of random points are generated. 
The detected galaxies are re-gridded into a rectangular box, 
the FKP estimate for the weighted density field is
\begin{eqnarray}
F(\bm{r}) = \frac{w(\bm{r})}{N}[n_{g}(\bm{r}) - \alpha n_{s}(\bm{r})],
\label{F}
\end{eqnarray}
where $w(\bm{r})=1/(1+\bar{n}(\bm{r})P)$ is the FKP weight, $\bar{n}(\bm{r})$ is the selection function (i.e. mean density) at the position $\bm{r}$, and P is an prior estimate of 
the power $P(k_1)$ at the scale of interest $k_1$, $n_{g}(\bm{r})$ and $n_{s}(\bm{r})$ refer to the number density 
of the observed galaxy catalog and the random mock catalog respectively, $\alpha$ is the real-to-mock ratio:
$\alpha = \sum w(\bm{r}_{\rm real})/\sum w(\bm{r}_{\rm random})$.
The  normalization factor $N$ in Eq. (\ref{F})  is given by
\begin{eqnarray}
N^2 = \int d^3r\bar{n}^2(\bm{r})w^2(\bm{r})
= \beta\sum_{\rm random}\bar{n}(\bm{r}_{i})w^2(\bm{r}_{i}),
\label{normal}
\end{eqnarray}
where  $\beta$ is the unweighted ratio of number of galaxies in the real (in our case simulation) catalog to that in the random mock catalog. 
In order to reduce the shot noise, the mock catalog is always set to contain much more galaxies than the real catalog. 
Here we choose $\beta$ to be 0.02. The power spectrum can then be estimated from the Fourier Transform of the weighted over-density field
$\langle |F(\bm{k})|^2\rangle$. The measurement error on power spectrum are estimated with $\Delta k/k = 0.125 $, for which the 
error of different $k$-bins can be regarded as uncorrelated \citep{1998ApJ...495...29G,1999MNRAS.304..851M}.

Our simulation data cube is in a pencil-beam shape, the amplitude and shape of the power spectrum estimated from the
cubic grid are biased, $\langle |F(\bm{k})|^2\rangle = W(k) \hat{P}(k)$. To correct for this effect, we compute the window function for this survey geometry by  
producing two  set of random catalogs, one is distributed only in the pencil-beam region, the other one 
in a cube region which encloses the pencil-beam region, the window function is then the ratio of the two power spectrum 
$P_{\rm pencil}(k)/P_{\rm cube}(k)$. We produce 10 pairs of samples and use the mean value to make the estimate. The true power spectrum is then obtained 
by dividing the window function. In summary, 
the power spectrum is obtained with the following steps:
\begin{enumerate}
\item Compute the selection function $\bar{n}(\bm{r})$ in each redshift bin. 
\item Produce the weighted over-density field $F(\bm{r})$) in the gridded box with grid spacing $1.0 \ensuremath{\,{\rm Mpc}}/h$ using the Nearest Grid Point assignment technique. 
\item Fourier transform the weighted over-density field, and compute the power spectrum $\langle |F(\bm{k})|^2\rangle$ with $\Delta k/k$ = 0.125. 
\item Correct the shape effect using the final window function,   $\hat{P}(k)= \langle |F(\bm{k})|^2\rangle /W(k)$. 
\end{enumerate}

The L-band receiver can cover a redshift range of up to $z=0.35$. In Fig. \ref{projected_error_galaxy_count} we show 
the projected error on power spectrum from a 20000-$\deg^2$ galaxy surveys for redshift $0.05<z<0.15, 0.15<z<0.25, 0.25<z<0.35$,
with integration time 48s, 96s, 192s and 384s per beam. On larger scales, the measurement precision is limited by 
the available number of modes (cosmic variance), while on the smaller scale it is limited by the available number of galaxies per 
cell (shot noise).  The best relative error in the power spectrum is achieved somewhere at $10^{-2} h/\ensuremath{\,{\rm Mpc}} <k < 10^{0} h/\ensuremath{\,{\rm Mpc}} $. 
At higher redshift the optimal point shift towards larger scales (smaller $k$), as the probed volume increases and 
the observed galaxy number density decreases.  At the BAO scale k $\approx$ 0.07 $h/\ensuremath{\,{\rm Mpc}}$, the signal to noise ratio can 
reach 5.0 at $z \approx 0.2, 0.25, 0.3$ and 0.35 respectively.   so in our Fisher estimation we will use these redshift regions 
as the survey volume of the FAST HI galaxy survey.  As can be seen from the figure, at the lower redshift, and especially on small $k$ values
the difference in the integration time per beam does not make much difference, however for the high redshift and large $k$ values the increase of 
integration time can significantly improve the measurement precision.  

The power spectra from the simulated observations are shown in Fig. \ref{fkp_multitime}. 
From top to bottom panel are the result with 48s, 96s, 192s and 384s respectively, each in the redshift range of 0.05-0.15 (left), 
0.15-0.25(middle), 0.25-0.35(right). In each subfigure, the theoretically expected power spectrum, the mock observation power spectrum, 
and the projected measurement error are plotted. For the noise, we also show separately the sample variance (corresponding to the size of the simulation
box) and the shot noise components. In this case, the cosmic variance is much larger as the size of simulation box 
spans  only a region of $10\times10\deg^2$.  

\begin{figure*}
  \centering
  \includegraphics[width=0.75\textwidth]{fkp_48s.pdf}
  \includegraphics[width=0.75\textwidth]{fkp_96s.pdf}
  \includegraphics[width=0.75\textwidth]{fkp_192s.pdf}
  \includegraphics[width=0.75\textwidth]{fkp_384s.pdf}
  \caption{The HI galaxy power spectrum from simulation. We show the theoretical matter power spectrum 
  (black dashed line),  the simulated measurement (blue solid line), the total projected error estimated with Fisher matrix 
  (red dash-dotted line) and the error due to sample variance (orange dash-dot-dotted line) and shot noise (green dotted line). }
  \label{fkp_multitime}
\end{figure*}


In this simulated galaxy survey, in almost all cases the measured power spectrum is larger than the theoretical power spectrum, 
as can be seen from Fig. \ref{fkp_multitime}. However, it does agree well with the projection. 
On small scales the shot noise dominated, which results in a nearly flat spectrum. This is particularly obvious at higher redshifts 
and for shorter integration times, where the number density of the detected galaxies is too small.  On larger scales, at lower redshifts 
and longer integration times, the shape of the power spectrum is more similar to the theoretical power spectrum, but there is still 
significant difference, and there are some large fluctuation at the large scales, thanks to the contribution of the cosmic variance. 
The overall normalization of the power is higher than the matter power spectrum (marked as theoretical), due to the fact that 
only the rare massive galaxies can be detected, and a clustering bias is introduced. The projected errors get larger in higher redshifts. 
However, we note that for a real survey with large sky area, the cosmic variance can be significantly reduced.

\subsection{Intensity Map Power Spectrum}

\begin{figure*}
  \centering
  \includegraphics[width=1.0\textwidth]{projected_error_im.eps}
  \caption{The projected error on power spectrum from FAST 20000 $\deg^2$ intensity mapping experiment with integration time of 48s, 96s, 192s and 384s per beam. At k $\approx$ 0.07 $h/\ensuremath{\,{\rm Mpc}}$ the S/N can reach 5.0 until redshift of 0.35, 0.55, 0.75 and 1.05 respectively.}
  \label{projected_error_im}
\end{figure*}

For intensity mapping survey, the measurement error of the power spectrum
can be written as \citep{2010ApJ...721..164S}

\begin{eqnarray}
    \frac{\sigma_{P}}{P}  & =  & 2\pi\sqrt{\frac{1}{V_{\rm eff}(k)k^2\Delta k}}, 
    \label{projectederror_im}
\end{eqnarray}
with the $V_{\rm eff}$ given in this case by
\begin{eqnarray}
V_{\rm eff}(\vec{k}) = V_{\rm sur}\left(1+\frac{\sigma_{\rm pix}^2V_{\rm pix}}{[\bar{S}(z)]W(\vec{k})^2P}+\frac{1/\bar{n}}{P}\right)^{-2}.
\label{Veff_im}
\end{eqnarray}
where $\bar{S}(z)$ is the average 21-cm emission flux density, $V_{\rm pix}$ is the pixel volume. 
The first term is due to sample variance, the second term is induced by the system thermal noise, 
and the last term is 
the shot noise due to the discreteness of the HI sources, with $1/\bar{n} \approx$ 100$h^{-3}\ensuremath{\,{\rm Mpc}}^3$ (\citet{2010ApJ...721..164S}). 
We model the angular resolution (the frequency resolution is much higher) as:
\begin{eqnarray}
W(k) = \exp\bigg[-\frac{1}{2}k^2r(z)^2\bigg(\frac{\theta_{\rm pix}(z)}{2\sqrt{2\ln{2}}}\bigg)^2\bigg]
\label{window}
\end{eqnarray}


\begin{figure*}
  \includegraphics[width=0.83\textwidth]{intensity_mapping_48s.eps}
  \vspace{-0.5ex}
  \includegraphics[width=0.83\textwidth]{intensity_mapping_96s.eps}
  \vspace{-0.5ex}
  \includegraphics[width=0.83\textwidth]{intensity_mapping_192s.eps}
  \vspace{-0.5ex}
  \includegraphics[width=0.83\textwidth]{intensity_mapping_384s.eps}
  \caption{The intensity mapping power spectrum for different integration time per beam from simulation.}
  \label{intensity_mapping_multitime}
\end{figure*}

The image cube is gridded into a rectangular box with grid spacing of 2 $\ensuremath{\,{\rm Mpc}}/h$,  here we take $2 \ensuremath{\,{\rm Mpc}}/h$ sized pixels as the 
standard sampling size, it is well below the BAO scale which is about $150 \ensuremath{\,{\rm Mpc}}/h$ and would not affect the cosmological result.  
The HI emission flux density field is related to the over density field by:
\begin{eqnarray}
\delta(\bm{r}) = S_{\rm HI}(\bm{r})/\bar{S}_{\rm HI}- 1,
\label{flux_to_density}
\end{eqnarray}
where $S_{\rm HI}(\bm{r})$ is the HI emission flux density at position $\bm{r}$ and $\bar{S}_{\rm HI}$ is the mean flux density of HI emission. 
We also add a Gaussian thermal noise of  0.74mJy, 0.52mJy, 0.37mJy and 0.26mJy per beam, corresponding to 48s, 96s, 192s and 384s 
integration time per beam. The simulated power spectrum is measured as follows:
\begin{enumerate}
\item Generating the mock data filled with Gaussian thermal noise,
\item Re-gridding the image data cube to rectangular box with grid spacing of 2 $\ensuremath{\,{\rm Mpc}}/h$,
\item Converting the flux density field to HI mass overdensity field, 
and Fourier transform the field to get the power spectrum $\langle |F(\bm{k})|^2\rangle$ with $\Delta k/k$ = 0.125.
\item Using the pencil beam survey window function to correct the shape effect.
\end{enumerate}

The statistical error is estimated as follows:
\begin{enumerate}
    \item Generating a synthetic product with a spatial resolution same as FAST beam size and a frequency resolution of 1MHz, S($r_{\rm x}$,$r_{\rm y}$,z),
    \item Generating a mock data filled with Gaussian thermal noise, N($r_{\rm x}$,$r_{\rm y}$,z), the noise data cube have the same size and 
    resolution as S($r_{\rm x}$,$r_{\rm y}$,z),
    \item Re-gridding the synthetic data S($r_{\rm x}$,$r_{\rm y}$,z) and the noise data N($r_{\rm x}$,$r_{\rm y}$,z) to a rectangular box with grid spacing of 2 $\ensuremath{\,{\rm Mpc}}/h$, obtaining S($r_{\rm x}^{\prime}$,$r_{\rm y}^{\prime}$,z) and N($r_{\rm x}^{\prime}$,$r_{\rm y}^{\prime}$,z), 
    \item Computing the mean signal of S($r_{\rm x}^{\prime}$,$r_{\rm y}^{\prime}$,z) and the variance of N($r_{\rm x}^{\prime}$,$r_{\rm y}^{\prime}$,z) at different redshift, obtain $\bar{S}(z)$ and $\sigma_{\rm pix}$, 
    \item Obtaining the projected error by use of Eq. (\ref{projectederror_im}).
\end{enumerate}

Figure. \ref{projected_error_im} shows the statistical error of the power spectrum that can be achieved by $20000 \deg^2$ FAST intensity 
mapping survey with integration time 48s, 96s, 192s and 384s per beam respectively. Here, we ignore for the moment the frequency range
of the L-band receiver system, but simply plot up to much higher redshifts, which may be accomplished with the wide-band single feed receiver, or a 
future UHF PAF receiver introduced in Section 2.1. For simplicity and easier comparison, here we also plot the projected error for the same integration 
time, though to the size of the beam and the number of available beams the required total survey time would be very different at the higher redshifts.

A comparison with Fig.~\ref{projected_error_galaxy_count} shows that for the lowest redshift bin ($0.05<z<0.15$), the projected errors 
are almost the same as the HI galaxy survey of the same integration time. We saw from Fig.~\ref{number_density_multitime} that even in this 
low redshift range, the comoving number density of the detected HI galaxies are decreasing with increasing redshift, showing that the FAST does not
detect all HI galaxies. Nevertheless, the detected HI galaxy has sufficiently high number density that it gives a good representation of the underlying 
HI galaxy sample and total mass density. However, for all the other redshift bins the projected error of the 
intensity mapping is much smaller than the HI galaxy survey. In fact, as redshift increases, the signal-to-noise ratio 
actually improves for a while as the survey volume increases and sample variance decreases, though 
eventually it begins to drop as the thermal noise becomes higher and HI signal becomes weaker.  
The shift in the optimal comoving scale (the minimum of the relative error in the power spectrum) at different redshifts are also 
much less than the HI galaxy surveys. 
At the BAO scale $k \approx$ 0.07 $h/\ensuremath{\,{\rm Mpc}}$, the S/N can reach 5.0 until
redshift of 0.35, 0.55, 0.75 and 1.05 for survey with 48s, 96s, 192s and 384s respectively. 

The simulated power spectrum measurement is shown in Fig. \ref{intensity_mapping_multitime}. The left three columns show the simulation with 
the same L-band receiver, while the column on the right shows a higher redshift ($0.35<z<0.45$) with the wide-band receiver system of a single feed. 
Compared with Fig.~\ref{fkp_multitime}, the thermal noise is replaced by the shot noise at the smaller scales.  Within the L-band (up to $z=0.35$), 
it is sub-dominant even for the shortest integration time considered here (48 s). So the intensity mapping can yield much nicer results than galaxy surveys on 
the small scales, and good precision can be achieved at  higher redshifts. On the larger scales where the cosmic variance dominates, 
the precision of this simulated measurement is limited by the size of the simulation box so we see some large fluctuations, 
but this can be significantly reduced with larger volumes.  The intensity mapping survey can efficiently map the  large-scale structure of 
BAO scale until $z \approx 0.35$ even with moderate integration time. Note that the error from intensity mapping is 
smaller at most of scales, but at small scales ($k \approx 1.0 h/\ensuremath{\,{\rm Mpc}}$) the noise is larger than HI signal at all redshifts due to the limited 
resolution.

For redshift $0.35<z<0.45$ case, i.e. the right column in Fig.~\ref{intensity_mapping_multitime},  the wide-band receiver is assumed to be 
used as the corresponding frequency range is beyond that of the L-band receiver.  
This receiver has a much higher noise level, the system temperature is 60 K, so the survey precision for the 
same integration time degraded a lot. Also, as it is a single 
feed receiver, the required survey time for accomplishing the same integration time would be much longer. In order to make a good intensity mapping survey 
at the higher redshifts, one must be able to have multi-beam receivers, preferably with low system noise. In the next subsection where the cosmological 
measurement is discussed, we shall assume that such a receiver, namely the UHF PAF receiver introduced in Sec.2.1 will be available in the future for 
such surveys. 



\subsection{Cosmological Constraints}

Because the limit of our survey volume, we estimate the induced measurement error on cosmological parameters by using the
Fisher matrix formalism \citep{1997PhRvL..79.3806T,2003ApJ...598..720S,2013MNRAS.430..747M}. The Fisher matrix for 
 parameter set $\{p_i\}$ is given by
\begin{equation}
F_{ij} = \int_{\vec{k}_{\rm min}}^{\vec{k}_{\rm max}}\frac{\partial \ln P(\vec{k})}{\partial p_{i}}\frac{\partial \ln P(\vec{k})}{\partial p_{j}}V_{\rm eff}(\vec{k})\frac{d\vec{k}}{2(2\pi)^3} 
\label{fisher}
\end{equation}
The usable range of $k_{\rm min}$ and $k_{\rm max}$ are assumed to be 
$k_{\rm min} = 10^{-3} h/\ensuremath{\,{\rm Mpc}}$ (from survey volume) 
and $k_{\rm max} = 0.1 h/\ensuremath{\,{\rm Mpc}}$. The observed power spectrum $P(\vec{k})$ is given by 
\begin{align}
    \begin{split}
    P_{\rm obs}(k_{\rm ref\perp},k_{\rm ref\parallel}) &= \frac{D_{A}(z)_{\rm ref}^2H(z)}{D_{A}(z)^2H(z)_{\rm ref}}b_{\rm HI}^2\left(1+\beta\frac{k_{\parallel}^2}{k_{\perp}^2+k_{\parallel}^2}\right)^2\\
                                                   &\hspace{1em}\times \left(\frac{G(z)}{G(z=0)}\right)^2P_{\rm m,z=0}(k) + P_{\rm shot},
    \end{split}
    \label{power_obs}
\end{align}
where $b_{\rm HI}$ is the linear bias factor of HI gas and the redshift space distortion factor $\beta = \Omega_{m}(z)^{0.6}/b_{\rm HI}(z)$. The Hubble parameter and the angular diameter distance can be computed for a model 
with dark energy equation of state parameterized in the form $ w(z) = w_{0}+w_{a}\frac{z}{1+z} $ \citep{2001IJMPD..10..213C}. 
To obtain useful constraints on cosmological parameters, it is necessary to break the degeneracy by combining the BAO data with data obtained from some other cosmological observations, e.g., CMB. The total Fisher matrix on distance parameters is given by 
\begin{equation}
F^{\rm tot} = F^{\rm CMB} + \sum_{i}F^{\rm LSS}(z_{i}), 
\label{fisher_tot}
\end{equation}
and $F^{LSS}(z_{i})$ is the Fisher matrix derived from the i-th redshift bin of the large scale structure survey. For FAST HI galaxy survey and intensity mapping observations, we divide the redshift region into several bins with equal redshift interval, which we set as 0.05. 

The HI gas mostly distributed in galaxies hosted by halos after the reionization, thus the HI bias can be modeled as the halo bias weighted by HI mass hosted by these halos:
\begin{align}
    b_{\rm HI} = \frac{1}{\rho_{\rm HI}(z)}\int_{M_{\rm min}}^{M_{\rm max}}dM\frac{dn}{dM}(M,z)M_{\rm HI}(M)b(M,z),
    \label{bias_HI}
\end{align}
where $dn/dM$ is the halo mass function, for which we use the fittings of (\citealt{2008ApJ...688..709T,2010ApJ...724..878T})
of the Sheth-Tormen function \citep{2002MNRAS.329...61S}.
The HI mass in a halo of mass $M$ is given by $ M_{\rm HI}(M) = A M^{\alpha}$ 
where the prefactor  $A$ will be canceled in the normalization of $\rho_{\rm HI}(z)$, 
and  $\alpha \simeq 0.6$ \citep{2015aska.confE..19S}. 
The halo bias in this ellipsoidal collapse model is modeled as,
$$
    b(M,z) = \frac{1}{\delta_{c}(z)}\bigg[\nu^{\prime2}+b\nu^{\prime2(1-c)}-\frac{\nu^{\prime2c}\sqrt{a}}{\nu^{\prime2c}+b(1-c)(1-c/2)}\bigg],
$$
where $a = 0.707$, $b = 0.5$, $c=0.6$, $\nu^{\prime} = \sqrt{a}\nu$, and  $\nu = \frac{1.686}{D(z)\sigma_{R}}$, 
where $D(z)$ is the linear growth factor \citep{2002MNRAS.336..112M}.

The bias for the HI galaxy survey and the HI intensity mapping survey is different. For the HI galaxy survey, the $M_{\rm min}$ is 
given by the minimum mass of detected galaxies in the survey, while for HI intensity mapping survey it is the minimum mass, about $\sim 10^6 \ensuremath{\, {\rm M}_{\odot}}$.





\subsection{Results}

We consider the constraint on dark energy EOS parameters ($w_{a}$,$w_{0}$) from the HI galaxy surveys and intensity mapping surveys. 
The cosmological constant point $w_0, w_a=(-1,0) $ is taken as the fiducial model. 
The L-band 19 beam receiver is assumed to be used for the redshift $0.05<z<0.35$. We then consider the cosmological constraint derived from (i) the 
L-band survey only; (ii) the L-band survey plus the existing wide-band single feed receiver; or (iii) the L-band plus the future
PAF receiver with 81-beams in the UHF band for higher redshift intensity mapping (see Sec.\ref{sec:instr} for descriptions). In each case, we consider 
the integration time per beam on the L-band is 48s, 96s, 192s, and 384s, respectively. As discussed in Sec.\ref{sec:instr}, the 48s integration can be 
completed in one scan of the 19-beam receiver.
When the wide band single feed receiver is used, for simplicity, we assume it acquires the same amount of integration time per beam. 
For the PAF in UHF band, we have considered instead two integration time per beam, 216 s corresponding to (one scan) or 432 s (two scans) of the 
PAF receiver. These are added to the L-band 192 s and 384 s respectively for illustration. 

These survey configurations are listed in the first column of Table \ref{covariance}. The precision of the dark energy EOS parameters $(\sigma_{w_0}, \sigma_{w_a})$ are given in the second column for galaxy survey and third column for 
intensity mapping survey. The required total observation time corresponding to each survey are 
given in the fourth column of the table, 
where the time required for each band are given separately.  
Note that the time listed is for observations, not counting offline time required for 
calibration, maintenance, etc., so the real time required to complete the survey would be even longer.  

We also show in Figure. \ref{w0_fld_wa_fld}  the error ellipses of the dark energy EOS parameters with the L-band 19 beam receiver and the wide band single 
feed receiver. The results for the L-band + PAF surveys are not shown as their error ellipses are much smaller.  
The DETF figure of merit, which is defined as the inverse of the area of the $2\sigma$ error ellipse \citep{2006astro.ph..9591A}, 
is shown in Figure. $\ref{w0_fld_wa_fld_FOM}$ for the L-band + wide-band receiver as well as the L-band + UHF-band PAF. 

\begin{table}
 \caption{Constraint on dark energy EOS parameters $(w_0, w_a)$ for FAST galaxy surveys (GS) and intensity mapping (IM) surveys  
 with a total of  $20000 \deg^2$ area and Planck prior. In first column the L denotes L-band 19 feed receiver, 
 w denotes wide band single feed receiver, and P denotes the PAF receiver in UHF band, and the number denotes 
 integration time in seconds per beam.  The last column shows the total time for the survey to be completed. }
 \label{covariance}
 \begin{tabular}{L{2cm}|C{1.4cm}|C{1.4cm}|L{2.3cm}}
 \hline
 Survey &GS $(\sigma_{w_0}, \sigma_{w_a})$ & IM $(\sigma_{w_0}, \sigma_{w_a)}$ & Observation Time (day)\\
 \hline
  L 48s & (0.46, 1.44)& (0.19, 0.53) &220\\
  L 96s & (0.33, 1.00)& (0.15, 0.43) &440\\
  L 192s & (0.25, 0.77)&(0.13, 0.36) &880\\
  L 384s & (0.17, 0.49)& (0.12, 0.33) &1760\\
 (L + w) 48s & (0.46, 1.44)& (0.18, 0.50) &220(L) + 2422(w)\\
 (L + w) 96s & (0.33, 1.00)& (0.14, 0.39) &440(L) + 4844(w)\\
 (L + w) 192s & (0.25, 0.77)& (0.11, 0.30) &880(L) + 9688(w)\\
 (L + w) 384s & (0.17, 0.49)& (0.09, 0.23) &1760(L) + 19376(w)\\
 L(192s) +P(216s) & - & (0.05, 0.12) &880(L) + 135(P)\\
 L(384s) +P(432s) &- & (0.04, 0.10) &1760(L) + 270(P)\\
 \hline
 \end{tabular}
\end{table}

From these we see the intensity mapping can achieve much higher precision in the measurement of the dark energy EOS parameters than the 
HI galaxy surveys. The intensity mapping survey with the shortest integration per beam (48s) has a figure of merit comparable with the 
HI galaxy survey of the longest integration time per beam (384s), but requires only 1/8 of the total observation time in the L-band. However, even with 
the intensity mapping survey, the figure of merit is only of order $10^1$, much less than the current optical surveys. This is not surprising as the L-band is
limited to relatively low redshifts ($z<0.35$). However,  we note that so far there has not yet been an HI survey providing constraints on the dark energy 
parameters, and the HI survey complements the optical survey as it uses a difference tracer, which would be valuable to reduce any possible 
systematic errors in the BAO measurement. 


\begin{figure}
  \centering
  \includegraphics[width=0.45\textwidth]{w0_fld_wa_fld.eps}
  \caption{The constrains on dark energy EOS parameters from FAST galaxies survey and intensity mapping, both combined with Planck CMB observation. The three colors are iso-probability contours for 0.68, 0.95 and 0.99, respectively. Galaxies survey (GS) is labeled as filled ellipse and Intensity mapping (IM) is labeled with solid lines.}
  \label{w0_fld_wa_fld}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.4\textwidth]{w0_fld_wa_fld_FOM.eps}
  \caption{The figure of merit for dark energy EOS parameters from FAST galaxies survey and intensity mapping combined with Planck CMB observation. Galaxy surveys (GS) are labeled as filled blue circles and intensity mapping surveys (IM) are labeled with filled red down-pointing triangle, for confidence limit of 0.95. The IM surveys with low frequency PAF are shown in green up-pointing triangle symbols.}
  \label{w0_fld_wa_fld_FOM}
\end{figure}

From Table \ref{covariance}, we also see that the constraints on the EOS parameters are only slightly improved by adding the surveys of equal integration
time per beam with the wide-band receiver. This is because this receiver has a higher system temperature (60 K) 
over its very wide frequency coverage. Furthermore, because it has only a single feed, to achieve such integration time per beam would require very long
observation time which are quite impractical. However, if equipped with a powerful PAF receiver, the measurement can be taken to much higher redshift 
in reasonable time--indeed, for the PAF parameter we assumed, it would take even less time to complete than the L-band. The figure of merit could then be 
lifted substantially, up to a level comparable with Dark Energy Task Force (DETF) stage IV experiments \citep{2006astro.ph..9591A}.  This shows that a
PAF receiver at the UHF band would be a very valuable addition to the FAST telescope. 









\section{Discussions}
\subsection{The choice of galaxy survey vs. intensity mapping}
Using the Fisher information, \citet{2018arXiv180906384C} develop a formalism to quantify the performance of galaxy redshift survey 
and intensity mapping when measuring large-scale structures. Under the assumption that the galaxy population follows 
a Schechter function form, 
\begin{eqnarray}
n(L) &=& \phi_* \left(\frac{l}{l_*}\right)^{\alpha} e^{-l/l_*},
\end{eqnarray}
the optimal strategy for survey can be found using the relative value of three parameters: 
$\{L_{\rm SN}$, $\sigma_{\rm L}$, $l_{\ast}\}$, where $L_{\rm SN}$ is the luminosity scale on which the voxels 
are susceptible to shot noise, $\sigma_{\rm L}$ refers to the rms noise per voxel, and $l_{\ast}$ is the characteristic luminosity. 
The $L_{\rm SN}$ is derived with $\sigma_{\rm SN}(L_{\rm SN}/l_{\ast}) = L_{\rm SN}/l_{\ast}$, where
\begin{eqnarray}
\sigma^{2}_{\rm SN}(l) = V_{\rm vox}\phi_{\ast}\int_{0}^{l}dl^{\prime}l^{\prime \alpha+2}e^{-l^{\prime}}.
\label{sigma_sn}
\end{eqnarray}
$V_{\rm vox}$ is the comoving volume of a voxel, $\alpha$ is the faint-end slope parameter of the luminosity function. 
Observations can be divided into  four limiting regimes for optimal strategy, as shown in Table (\ref{regime}). 
In regime 1, the instrument noise is much smaller than $l_{\ast}$, and confusion effect is small, galaxy detection is optimal. 
In regime 2, the optimal strategy is somewhere intermediate between the intensity mapping and galaxy detection, 
because the voxels with $L \geqslant \sigma_{\rm L}$ will suffer from confusion noise. In regime 3, 
the instrument noise in a voxel is very large, $l_{\star} < \sigma_{\rm L}$, the intensity mapping will be the only choice. 
Regime 4 corresponds to a large effective number of galaxies per voxel, galaxy detection will suffer from large confusion noise, intensity mapping is optimal. 

For the FAST HI survey, assuming a voxel of an angular resolution of $0.08^\circ$ and bandwidth of 200km/s, we compute the 
evolution of $L_{\rm SN}/l_{\ast}$ and $\sigma_{\rm L}/l_{\ast}$ in redshift and show the results in 
Figure. \ref{lsn_sigmal_l}. We use the luminosity function given by \citet{2003ApJ...592..819B}, 
with the following parameters: 
$\phi^{\ast} = 5.11\times10^{-3}h_{70}^3$ Mpc$^{-3}$, $\log (L_{\ast}/L_{\odot})=10.36+\log h_{70}$ and $\alpha= -1.05$. 
The redshift evolution of Schechter luminosity function parameters is modeled as\citep{1999ApJ...518..533L},
\begin{eqnarray}
     \alpha(z) &=& \alpha(z_{0}), \label{lf_evolution_alpha}\\
     M^{\star}(z) &=& M^{\star}(z_{0}) - Q (z - z_{0}), \label{lf_evolution_M}\\
     \varphi^{\star}(z) &=& \varphi^{\star}(z_{0})10^{0.4P (z - z_{0})}, \label{lf_evolution_phi}
\end{eqnarray}
where $P = 1.0$ and $Q = 1.03$ \citep{2015MNRAS.451.1540L} in the $r$-band. 

\begin{table}
 \caption{Four limiting regime defined by the relative value of the luminosity scale where the voxels are
  highly susceptible to shot noise, $L_{\rm SN}$, the rms noise per voxel, $\sigma_{\rm L}$ and the characteristic
   luminosity for a certain luminosity function, $l_{\star}\}$.}
 \label{regime}
 \begin{tabular}{ccc}
  \hline
  number&regime&optimal strategy\\
  \hline
  1&$L_{\rm SN} < \sigma_{\rm L} < l_{\star}$ & galaxy detection\\
  2&$\sigma_{\rm L} < L_{\rm SN} < l_{\star}$ & galaxy detection/intensity mapping$^a$\\
  3&$L_{\rm SN} < l_{\star} < \sigma_{\rm L}$ & intensity mapping\\
  4&$l_{\star} < L_{\rm SN}$ & intensity mapping\\
  \hline
  \multicolumn{3}{l}
  {\parbox{8cm}{$^a$ Here the optimal strategy is an intermediate between the  intensity mapping and galaxy detection observables.}}\\
 \end{tabular}
\end{table}

Figure. \ref{lsn_sigmal_l} shows the redshift evolution of $L_{\rm SN}/l_{\ast}$  (top panel), and $\sigma_{\rm L}/l_{\ast}$ (bottom panel).
We see $L_{\rm SN}\leqslant l_{\ast}$ at redshift $\lesssim$ 0.8. and $\sigma_{\rm L}\leqslant l_{\star}$ 
at redshift $\lesssim$ 0.13, 0.16, 0.19 and 0.23 for integration time of 48s, 96s, 192s and 384s respectively. 
Combining with Table (\ref{regime}), it shows for FAST the galaxy redshift survey is the optimal strategy at 
redshift $\lesssim$ 0.13, 0.16, 0.19 and 0.23 for survey with 48s, 96s, 192s and 384s integration time per beam. 
This is in agreement with what the projected error shown in Figure. \ref{projected_error_galaxy_count} and 
Figure. \ref{projected_error_im}.  We note that the voxel size in the projected error calculation for the 
intensity mapping is 2 $\ensuremath{\,{\rm Mpc}}/h$, which is larger than the voxel size for galaxy detection. This makes the redshift points 
 where intensity mapping is better than galaxy detection a little lower than the redshifts shown in Figure. \ref{lsn_sigmal_l}.


\begin{figure}
  \centering
  \includegraphics[width=0.38\textwidth]{Lsn.eps}
  \includegraphics[width=0.38\textwidth]{sigma_l_l.eps}
  \caption{The redshift evolution of $L_{\rm SN}/l_{\ast}$ and $\sigma_{\rm L}/l_{\star}$, assuming an angular resolution of $0.08^\circ$ 
   and a bandwidth of 200km/s. The jump at redshift 0.35 in the bottom panel is due to the change of receiver system.}
  \label{lsn_sigmal_l}
\end{figure}

\subsection{Foreground}
One of the most challenging problems in intensity mapping experiment may be the contamination from the foreground radiation, which is 
several orders of magnitude larger in amplitude than the HI intensity signal.  It can in principle be subtracted, and the true signal recovered, 
based on the fact that the frequency dependency and some statistical properties of the foreground are different form the true signal. 
Sophisticated mathematical methods have been developed  \citep{2011PhRvD..83j3006L,2013MNRAS.434L..46S,2013MNRAS.429..165C,2014MNRAS.441.3271W,2014ApJ...781...57S,2015MNRAS.447..400A,2016ApJS..222....3Z}.
The intensity mapping experiment with the Green Bank Telescope (GBT) (\citet{2013MNRAS.434L..46S}) have shown that the 
 foregrounds can indeed be suppressed significantly,  though at present a positive detection of 21cm auto-power spectrum is yet to be achieved. 
Here we assume that after  a successful foreground subtraction, the contamination can be reduced to 
the thermal noise level \citep{2015MNRAS.454.3240B}. 




To investigate the impact of foreground on the FAST IM survey, we made a simple 
test of foreground removal in our simulation. 
We produce the foreground with the global sky model (GSM) \citep{2008MNRAS.388..247D,2017MNRAS.464.3486Z}, 
convolved with a frequency-dependent beam and add to the noise-filled data cube. Then using a third order log-log
 polynomial fitting, we find the foreground can be removed effectively. The residue difference power spectra between 
 the original 21cm signal and the one obtained by removing the simulated foreground  with the polynomial fit 
 in the image cube are shown in Fig. \ref{spectra_with_foreground} with different integration times per beam for a Gaussian beam.  The power spectra show that the differences are small, at about 1 percent level at k $\approx$ 0.1 $h/\ensuremath{\,{\rm Mpc}}$ for these conditions. As expected, the amplitude of this error power spectra generally decreases  as integration time increases.
The differences at larger scales (smaller k) is larger, because a few large scale modes (modes with small k$_{\parallel}$) are removed by the foreground removal methods, 
in principle such large scale modes can be reconstructed via cosmic tidal reconstruction \citep{2018PhRvD..98d3511Z}. 
Using the foreground removed cube, we obtain the same constraints on the 
cosmological parameters in high accuracy. 


\begin{figure}
  \centering
 \includegraphics[width=0.45\textwidth]{ps_foreground_compare.eps}
   \caption{The relative error power spectra, defined by $(P_{\rm HI}-P_{\rm fr})/P_{\rm HI}$, where $P_{\rm HI}$ is the true 21cm power, 
  and $P_{\rm fr}$ is the foreground removed power  at different redshifts and integration time.}
  \label{spectra_with_foreground}
\end{figure}
\begin{figure}
  \centering
 \includegraphics[width=0.45\textwidth]{bessel_spectra.eps}
   \caption{Same as Fig.\ref{spectra_with_foreground}, but the beam is modeled by Bessel function.}
  \label{spectra_with_foreground_bessel}
\end{figure}


Besides the Gaussian beam profile, we also made a test for the case of oscillatory side lobes. 
 In Fig.\ref{spectra_with_foreground_bessel} we make the same foreground subtraction excise, but for a beam function modeled
 using the Bessel function as $|2\times J_1(x)/x|^2$, where $x=3.23\times\theta/\sigma_b$. The result is qualitatively similar to the Gaussian case,
 though there are some differences. 


The simulation of the foreground and its removal presented above may still
be too simplified. In reality, the the beam response
could be more complicated, and we only have an imperfect knowledge which must itself be determined from observation, 
and there are also irregularities in the beam and bandpass, polarization leakage, the $1/f$ noise and variations of the system gain, etc., making it much harder to remove. So the actual impact of the foreground could be higher than this simple estimate. 
 Obviously, the difficulty of the foreground removal depends on the design of the telescope. The stability, regularity and dynamic
range all affects the induced foreground contamination. A detailed study of foreground subtraction requires a realistic assessment of FAST telescope response,
which could only be obtained with actual observational data. This is beyond the scope of the present paper. However, the high sensitivity of the FAST is advantageous for 21cm extraction: the individual voxel signal-to-noise ratio is relatively high, so it 
is easier to be detected.

\section{Summary}
In this paper, we make a detailed study of large area drift scan HI survey with the FAST telescope. We considered using the existing L-band 19 beam 
receiver, the single feed wide band receiver for lower frequency (UHF), and also contemplated using a future UHF-band PAF receiver. 
We simulated observation of HI galaxies,  the number density of the detectable HI galaxies decreases rapidly as the redshift increases, and also 
due to the larger beam size and smaller galaxy size at higher redshifts, the mean number of galaxies within each voxel which we called confusion rate 
also increases. For the detected galaxies, there are also  measurement errors in both the direction of the line of sight and the direction perpendicular to it. 
We estimated such errors, but found that the main source of HI galaxy power spectrum error still comes from the shot noise on small scales and 
cosmic variance on large scales. We also considered the HI intensity mapping observation. 

The precision of the power spectrum measurement 
are forecasted using Fisher matrix for large survey areas, and we also make mock observations for both galaxy survey and intensity mapping survey
using simulation. With intensity mapping the power spectrum can be measured with high precision.  We find that the 
FAST can effectively detect the individual galaxy till z $\approx$  0.2, 0.25, 0.3 and 0.35 or map the large-scale structure 
with intensity mapping till z $\approx$ 0.35, 0.55, 0.75 and 1.05 respectively, if we assume 48s, 96s, 192s and 384s integration time per beam. 
Generally the HI intensity mapping observations can yield much more precise measurement, though 
the HI galaxy survey can also achieve nearly optimal measurement at lower redshifts,  with $z \lesssim 0.13, 0.16, 0.19$ and 0.23 for surveys 
with 48s, 96s, 192s and 384s integration time per L-band beam.  

We find that the FAST HI intensity mapping survey can produce a good measurement of the underlying power spectrum, and use the BAO method to 
measure the dark energy equation of state parameters. Such a measurement with a radio tracer is complementary to the optical BAO measurements, 
and reduce possible systematic errors. We also made a study of the impact of foregrounds on the measurement by simulation, and found that
it does not significantly affect the results, at least under the simplified conditions assumed in the study. 

With existing receivers the achievable precision is moderate, as it is largely limited to $z<0.35$, for higher redshifts 
the wide-band single feed would take too long time to complete the survey. However, if equipped with a UHF-band multi-beam receiver, higher redshifts 
can be observed more efficiently. We considered the case of a UHF-band PAF receiver with frequency coverage of $0.5 \sim 1.0\ensuremath{\, {\rm GHz}}$, and found that
the resulting survey may yield dark energy figure of merit up to $10^2$, comparable with the DETF stage IV results. 
This shows that a state-of-the-art PAF multi-beam receiver would be a very valuable addition to the FAST telescope.

\section*{Acknowledgements}
This work was first inspired by the late Prof. Rendong Nan. The parameters of the FAST receivers are provided to us by 
Di Li,  Youling Yue, Lei Qian and Jin Fan.  We also thank Shifan Zuo and Yichao Li for discussions.
This work is supported by the National Natural Science Foundation of China (NSFC) key
project grant 11633004, the Chinese Academy of Sciences (CAS) Frontier Science Key Project 
QYZDJ-SSW-SLH017 and the CAS Interdisciplinary Innovation Team grant (JCTD-2019-05), the NSFC-ISF joint research program No. 11761141012,  the MoST grant 2016YFE0100300, 
 the CAS Strategic Priority Research Program XDA15020200, the National Key R\&D Program 
2017YFA0402603, and the NSFC grant 11773034.


\bibliographystyle{mnras}
