https://orcid.org/0000-0001-6594-107X
ABSTRACT
Neutral hydrogen (H i) serves as a competitive tracer of the large scale structures, especially with the advent of more intensity mapping H i surveys. In this work, we present a simulation-based framework to forecast the H i power spectrum on non-linear scales (|$k\gtrsim 1\ {\rm Mpc^{-1}}$|), as measured by interferometer arrays like MeerKAT in the low-redshift (|$z\le 1.0$|) Universe. Building on a galaxy-based H i mock catalogue, we meticulously consider various factors, including the emission line profiles of H i discs and some observational settings, and explore their impacts on the H i power spectrum. We find that the H i power spectrum is relatively insensitive to the profile shape of H i emission line at these scales, while showing a strong correlation with the profile width. We propose an empirical model to simulate the emission line profile width for each H i source. The resulting H i power spectrum is consistent with the results from the IllustrisTNG hydrodynamical simulation and follows the trend of the measurements obtained by MeerKAT at |$z\approx 0.44$|, though with a significantly lower amplitude. We demonstrate how the H i abundance |$\Omega _{\rm HI}$| and the amplitude parameter in our width model can be constrained with the MeerKAT measurements, though a strong degeneracy is uncovered. Our work shows the potential to constrain statistical properties of H i emission line profiles with future H i intensity mapping experiments.
1 INTRODUCTION
In a low-redshift Universe (|$z\le 1.0$|), more than 95 per cent of atomic hydrogen (H i) gas resides in the cold dense region of galaxies within dark matter haloes, where it can be self-shielded from ultraviolet photons (see e.g. Krumholz, McKee & Tumlinson 2009; Diemer et al. 2018; Villaescusa-Navarro et al. 2018). The detectable 21 cm radiation emitted by the hyperfine transition of H i makes it a novel and competitive probe to measure the matter distributions and structures in the post-Epoch of Reionization (post-EoR) Universe at various spatial scales.
In the very local Universe (|$z\le 0.1$|), large dish radio telescopes such as Arecibo, Parkes, and Five-hundred-metre Aperture Spherical radio Telescope (FAST) have detected 21 cm signals emitted by extragalactic sources (Barnes et al. 2001; Haynes et al. 2018; Zhang et al. 2024). The extragalactic H i surveys performed by these telescopes have provided a wealth of combined information on H i mass, redshift, velocity profile, and yielded valuable insights into various aspects of H i statistics, including the H i mass function (HIMF) and the consequent H i abundance |$\Omega _{\rm HI}$| (e.g. Martin et al. 2010; Jones et al. 2018), the H i-halo mass relation (Guo et al. 2017, 2020), etc. However, at higher redshifts (|$z \ge 0.1$|), as the distances between the observer on Earth and the H i sources continually increase, the emission lines become significantly fainter and harder to distinguish from the noise. To overcome this limitation, the stacking method and intensity mapping (‘IM’ hereafter) technique is widely applied, both of which can measure the combined signal emitted by several unresolved H i sources instead of resolving individual objects, but only IM retains the spatial information (e.g. Battye, Davies & Weller 2004; Chang et al. 2008; Mao et al. 2008; Wyithe & Loeb 2009; Battye et al. 2013).
Numerous radio telescopes utilizing IM are currently either under construction or already in operation, such as the BAO from integrated neutral gas observations (BINGO; Wuensche and the BINGO Collaboration 2019), the Canadian Hydrogen Intensity Mapping Experiment (CHIME; Bandura et al. 2014), the Square Kilometre Array (SKA; Dewdney et al. 2009), and its pathfinder MeerKAT (Jonas & MeerKAT Team 2016). These telescopes fall into two main categories: single-dish telescope and interferometer array. In the case of the MeerKAT and SKA-mid arrays, they can be used in both observation modes, providing flexibility. The single-dish mode is optimized to observe comparatively large patches and collect the overall signal from relatively fast scans. On the contrary, the interferometer mode captures the signals from smaller fields of view by combining multiple small dishes. With the baselines, i.e. the distances between dishes ranging from metres to kilometres, the interferometer mode can achieve higher resolution, allowing for detailed studies of non-linear scales.
Notably, Paul et al. (2023) presented the first measurements of the H i power spectrum on non-linear scales obtained by the interferometer mode of MeerKAT, making it urgent to obtain a theoretical reference on these scales. Therefore, this work mainly focuses on the H i power spectrum at non-linear scales, to provide a finer theoretical framework for upcoming and existing H i IM surveys with interferometers and improve the interpretation of the data.
When dealing with non-linear scales, the power spectrum of compact astrophysical objects such as galaxies or H i sources is not solely influenced by the linear part. The shot noise term, which arises from the discrete characteristics of these compact objects, commonly behaves like a plateau and dominates as k increases to non-linear scales (e.g. see Wolz et al. 2019). Note that large scales in real space correspond to small k modes in Fourier space, and vice versa. However, in the case of H i, the significant internal movement of the H i gas within the discs, including both circular rotation and random motion, results in the elongation of the H i sources along the line of sight (LOS) due to the Doppler effect in redshift space. Therefore, since the H i discs are no longer compact (though only in the LOS), the H i power spectrum appears to continuously decrease with increasing k compared to a constant shot noise term (Sarkar & Bharadwaj 2019; Zhang et al. 2020; Padmanabhan et al. 2023).
Addressing the kinematics of the H i gas is crucial when measuring the power spectrum on small scales. Previous studies using the particle data in hydrodynamical simulations show no ‘shot noise plateau’ for H i power spectrum in redshift space (Villaescusa-Navarro et al. 2014, 2018). Other studies have constructed some models for the H i velocity profiles, confirming the suppression of power spectrum at small scales by the inner motion of the H i discs (see e.g. Sarkar & Bharadwaj 2019; Zhang et al. 2020). However, these models have yet to fully explore certain aspects, such as the choice of halo- or galaxy-based H i distributions, reasonable H i disc sizes, the impact of velocity profile shape on the H i power spectrum, and so on. The observational settings in the power spectrum measurements are also widely overlooked.
To address these questions and refine our theoretical framework to better align with measurements, we employ the neutraluniversemachine galaxy-H i catalogue, which is built on the top of universemachine simulation (Behroozi et al. 2019) and derived using the latest empirical model proposed by Guo et al. (2023). In the empirical model, H i statistics, e.g. H i mass function (HIMF), H i stellar mass relation, and H i abundance |$\Omega _{\rm HI}$|, are well constrained by existing observational measurements (see more details in Section 2).
Given that the catalogue lacks particle information, and therefore velocity information, it is critical to model and assign velocity profiles to each source. We explore the impact of its shape on the H i power spectrum by testing three toy models. Then we construct a model for the velocity dispersion, and use some statistics of emission line profile measured by extragalactic H i surveys to constrain the parameters in our model. The predicted 21 cm brightness power spectrum is calculated as equation (1) in this work.
where |$\overline{T}_b$| is the mean brightness temperature; |$P_{\rm HI}(k,\mu , z)$| is the cylindrical power spectrum of normalized H i density fluctuation |$\delta _{\rm HI}=\rho ({\bf x})/\overline{\rho }-1$|. |$\overline{T}_{\rm b}$| can be determined as follows (Berti, Spinelli & Viel 2023),
where |$\Omega _{\rm HI}$| is the H i abundance, |$H_0\equiv 100h\ {\rm km\ s^{-1}}{\rm {Mpc}^{-1}}$| is the Hubble constant, |$H(z)$| is the Hubble parameter. By normalizing the H i abundance in |$P_{\rm HI}$| and incorporating it into |$\overline{T}_{\rm b}$|, this work separately presents the H i bias |$b_{\rm HI}$| and H i abundance |$\Omega _{\rm HI}$|. Lastly, we illustrate how different observational settings affect the H i power spectrum and constrain our model based on the measurements.
This paper is organized as follows: Section 2 provides a brief introduction to the simulations and measurements used in this work. Section 3 describes how we simulate the emission line profile based on H i and subhalo/galaxy catalogues. Section 4 presents the results, including the cylindrical 2D and spherically averaged 1D power spectra in redshift space, H i bias, and a discussion on the shot noise term. Section 5 examines how observational settings impact the measured power spectrum and the constraint results. Finally, Section 6 summarizes the conclusions of this work. We adopt a flat, |$\Lambda$|CDM cosmology with parameters (|$\Omega _m=0.307, \Omega _{\Lambda }=0.693,h=0.678,\sigma _8 = 0.823, n_s = 0.96$|) consistent with Planck 15 results (Planck Collaboration XIII 2016).
2 OVERVIEW OF THE DATA
In this section, we give a brief introduction to the mock H i catalogues applied in this work including neutraluniversemachine and IllustrisTNG100-1, and also the measurements of the 21 cm power spectrum from the MeerKAT telescope. Most of our work is based on neutraluniversemachine catalogue, so all results are based on this, unless otherwise specified.
2.1 Modelling H i catalogues
2.1.1 neutraluniversemachine catalogue
neutraluniversemachine (hereafter ‘NUM’) is an empirical model framework proposed by Guo et al. (2023), which can simultaneously give self-consistent predictions for the evolution of H i and molecular hydrogen (H|$_2$|) in a redshift range of |$0\le z\le 6$|. To achieve this, they built the model as functions of redshift z and properties of (sub)haloes and galaxies (see their equations 1 and 15) and jointly constrain the 15 free model parameters using observational measurements in the redshift range of |$0\lt z\lt 6$|. The best-fitting model successfully describes the H i–halo relation, as well as the evolution of the cosmic H i abundance. Explicitly, the H i mass within a halo or subhalo has been found to be correlated with the virial halo mass |$M_{\rm vir}$| (Battye et al. 2013; Villaescusa-Navarro et al. 2018; Guo et al. 2020), the halo formation time |$z_{\rm form}$| (Guo et al. 2017; Li, Guo & Mao 2022), and the star formation rate (hereafter ‘SFR’) (Guo et al. 2021; Saintonge & Catinella 2022), which motivates the empirical functional form of the H i gas in the neutraluniversemachine model.
In this paper, we use the public H i catalogues of the neutraluniversemachine model, which is derived from the public catalogues of universemachine DR1. These catalogues are generated by running the empirical model universemachine (see Behroozi et al. 2019) in the Bolshoi-Planck N-body simulation (Klypin et al. 2016), with a box size of 250 Mpc h−1 and a dark mass particle resolution of |$2.3\times 10^8\ M_{\odot }$|.
2.1.2 IllustrisTNG100-1 catalogue
The IllustrisTNG (hereafter ‘TNG’) simulation (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018, 2019; Pillepich et al. 2018; Springel et al. 2018) is a set of state-of-the-art hydrodynamical simulations, containing a wide range of physical processes of galaxy formation. In this work, we adopt the simulation set of TNG100-1 with a box size of 110.7 Mpc and a dark matter particle mass resolution of |$7.5\times 10^6M_{\odot }$|.
To be more accurate, we employ the updated H i mass calculated by Diemer et al. (2018). In previous works, the ultraviolet photons emitted from quiescent gas cells with density below a so-called star-formation threshold are usually ignored (Marinacci et al. 2017) or set to the cosmic UV background (Lagos et al. 2015) while Diemer et al. (2018) consider the contribution from these cells to calculate the atomic-to-molecular transition in TNG simulation. There are several models discussed and updated in this paper, and for simplicity, we employ the ‘K13’ model, which is an analytical model to predict the fraction of molecular hydrogen H2 and H i on the top of neutral gas surface density and metallicity (Krumholz 2013).
2.2 Autocorrelation H i intensity mapping measurement
Paul et al. (2023) (hereafter P23) reported the first detection of the non-linear H i power spectrum independently at redshifts |$z=0.32$| and |$z=0.44$| using the L-band (900–1670 MHz) data obtained by the interferometer mode of MeerKAT radio telescope. As the precursor to SKA radio arrays, the MeerKAT radio telescope can be used in two different modes (single-dish mode and interferometer mode) for different cosmology purposes (Square Kilometre Array Cosmology Science Working Group 2020). MeerKAT consists of 64 |$13.5\,$|m dish antennas located in South Africa.1
The measurements of the H i power spectrum came from the Deep2 field (|$\alpha =04^{\rm {h}}13^{\rm {m}}26.4^{\rm {s}}, \delta =-80^{\circ }00^{\prime }00^{\prime \prime }$| in the Southern Hemisphere), in which only limited foreground radio sources exist. The visibilities were calibrated via the processMeerKAT (Collier et al. 2021) software pipeline. They chose two RFI-free frequency ranges centred at 986 and 1077.5 MHz with a bandwidth of |${\sim }\, 46$| MHz, equivalent to |$\Delta z \sim 0.06$|, small enough to ignore the cosmological evolution. The foreground avoidance technique was applied in the analysis, and the frequency resolution of the data is 209 kHz, equivalent to a comoving distance |$0.9596\ \mathrm{Mpc}$| at |$z=0.32$| and |$1.065\ \mathrm{Mpc}$| at |$z=0.44$| (see calculation details in Section 5.1.2).
3 SIMULATED EMISSION LINE PROFILE
H i emission line profiles have been widely detected by extragalactic H i surveys in the local Universe (|$z\le 0.1$|) (Barnes et al. 2001; Haynes et al. 2018; Maddox et al. 2021; Zhang et al. 2024). At higher redshifts, most H i sources are too distant to be identified by current radio telescopes, which falls short of the completeness requirements needed for common statistical studies.
When applying the intensity mapping technique, such selection effects are eliminated, and the velocity dispersions in the combined signals can still be resolved by some interferometer arrays like MeerKAT and therefore impact the measured H i power spectrum.
To simulate the emission line profile, the direct way is to run hydrodynamic simulations which include the velocity information of particle data (e.g. Villaescusa-Navarro et al. 2018; Davé et al. 2019). However, hydrodynamic simulations themselves are typically computationally demanding, and the particle H i data are usually obtained during the post-processing, making them heavily model-dependent and difficult to parametrize. Another method, which applies mock galaxy-based or halo-based H i catalogue, is more widely adopted, as it is in this work.
Using mock catalogues is more convenient and reasonably reliable, since in the low redshift Universe, almost all H i reside in galaxies and haloes (Villaescusa-Navarro et al. 2018). Based on different N-body dark matter simulations, both Sarkar & Bharadwaj (2018, 2019) assume that H i resides in haloes and has different velocity profiles as functions of the properties of their host dark matter.
Sarkar & Bharadwaj (2019) constructed a double-horn-like model for the velocity profile as a function of halo circular velocity |$v_{\rm cir}$|. They found that including the inner motion of the H i disc could considerably enhance the effect of the FoG, that is, suppress the H i power spectrum in large k modes, and that the suppression can be modelled by a Lorentzian damping profile. Zhang et al. (2020) adopted the relation between H i velocity dispersion and the host halo mass found by Villaescusa-Navarro et al. (2018) (see equation 26). They constructed a Gaussian model for the H i velocity profile as a function of the host halo mass.
Nevertheless, the plausibility of using halo-based models without considering subhalo/galaxy structures to predict H i power spectra needs to be deliberated, especially on non-linear scales, and the shape and width of H i emission line profile are modelled quite differently in these works. In this section, we verify that both galactic H i catalogues and an accurate velocity profile are needed to model the real signal and we give a prescription for how to incorporate the velocity profile with the calculation of the H i power spectrum.
3.1 Which elements matter?
To simulate the H i signal for a certain survey, several factors need to be carefully considered, which include the diameter of H i discs, stretching length of H i velocity dispersion in redshift space, the distances between galaxies. We can ask the following questions to elicit the key point. Can the resolution of the H i survey resolve the structures within haloes? Does the velocity dispersion caused by the Doppler effect have a non-negligible effect on the power spectrum? Is the H i density profile necessary to be modelled?
Taking the observation by P23 as an example, Fig. 1 transforms all the factors into length for a better comparison and articulates which elements matter for the observation. Upper panels depict lengths along the LOS direction, accounting for Redshift Space Distortion (RSD) effects, while the lower panels are lengths perpendicular to the LOS direction. The dashed grey lines in all panels are the length resolution, which is either the frequency resolutions (209 kHz) along LOS or the largest k modes in the measurements perpendicular to LOS (|$k \lesssim 10.5\ {\rm Mpc}^{-1}$| at |$z=0.32$|, |$k \lesssim 11.5\ {\rm Mpc}^{-1}$| at |$z=0.44$|).
The probability distribution of the scales for H i profile (black lines), H i emission line width in different direction (red lines), and the gap between galaxies (blue lines) in |$k_{||}$| direction (top panels) and in |$k_{\perp }$| direction (bottom panels). The grey dashed lines in the top panels represent the frequency resolution, i.e. channel width adopted by the autocorrelation measurement (209 kHz in |$z=0.32, 0.44$|). In the bottom panels, the grey dashed lines represent the scale corresponding to the largest k modes (|$k \lesssim 10.5\ {\rm Mpc}^{-1}$| at |$z=0.32$|, |$k \lesssim 11.5\ {\rm Mpc}^{-1}$| at |$z=0.44$|) shown in MeerKAT measurements. For greater visualization, lengths smaller than the frequency resolution or the selected k-mode are shaded in grey, where the length is less likely to be resolved by the data.
The right white region of the grey lines represents lengths that can be resolved by the data, while left shadow regions represents length smaller than the resolution. Imagine a mesh field gridded by the resolution lengths within which H i sources are placed. If sources are precisely located at the centre of each grid, then all length smaller than the resolution cannot be resolved. However, in reality, many sources are situated at random positions within the grids. Therefore, not all lengths smaller than the resolution are necessarily unresolved by the telescope, as not all sources are exactly at the centre of the grid. Therefore, lengths slightly smaller than the resolution are still likely to be resolved. It is tricky to specify a particular length value below which the factor and corresponding lengths can be neglected, so we recommend to disregard factors of lengths much smaller than the resolution and simulate the remaining factors accurately.
In this figure, the black lines illustrate the histogram of diameters of H i density profile that are estimated using the H i size-mass scaling relation obtained by Wang et al. (2016); red lines is the histogram of the simulated emission line width denoted as ‘|$W_{\rm 50}$|’ obtained by fitting to the width function of |$W_{\rm 50}$| measured by ALFALFA survey (see more details in Section 3.3.2); blue lines represent the histogram of ‘galaxy gap,’ which is the distances between galaxies inside the same host haloes in redshift or real space. Notably, the blue histograms differ considerably along the LOS and vertically because of the peculiar velocity of galaxies.
There are three points that can be confirmed through Fig. 1. First, it is necessary to adopt galaxy-based H i catalogues, because the telescope can resolve the galactic structures inside haloes. As shown in all of the panels, the majority of blue histograms are in the right-hand of the dashed lines. Secondly, it is unnecessary to consider the H i density profiles inside the H i discs. The sizes of H i discs are much smaller than the scale of our concerns in both directions, so it is safe to neglect the H i density profiles and to regard the H i discs as point masses in this work. Thirdly, although only a small part of the velocity dispersion along the LOS direction is larger than the frequency resolution, it is essential to account for this dispersion in the computation. This is because the larger velocity dispersion, the larger H i mass (e.g. see fig. 3 in Oman 2022) and thus including the velocity dispersion can severely suppress the shot noise term [see equation (11), where the large mass term accounts for a large fraction].
3.2 Calculation methods in simulations
We outline the steps we take to generate a H i density map in redshift space and compute the H i power spectrum as follows.
1. Calculate the H i position in redshifts space and set periodic boundary conditions;
where the |${\bf x}$| is the comoving position of H i sources in real space; |${\bf s}$| is the comoving position in redshift space; |$v_\parallel$| is the peculiar velocity projected in line of sight (LOS) direction. Plane-parallel approximations are adopted here; |$H(z)=H_0\cdot E(z)$| is the Hubble parameter.
2. Evenly cut the 2503 Mpc3 h−3 simulation box into 12003 grids. If frequency resolution is considered, the grid number along this axis should be reset accordingly. Then locate the centre of each H i source in the mesh according to their position |${\bf s}$|.
3. Assign a single-peak or double-horn emission line profile along the line of sight (LOS) for each H i source, with the profile width, such as |$W_{\rm 50}$| or |$W_{\rm 20}$| in extragalactic surveys, being solely determined. In this work, we use a Gaussian function for the shape and full-width at half-maxima (FWHM) for the width, also known as |$W_{\rm 50}$| (see more details in Section 3.3.2). The standard deviation |$\sigma$| of the Gaussian profile can be calculated accordingly.
Along the chosen LOS direction, select the grids within the |$\pm 2\sigma$| length near the ith source and assign its mass |$M_i$| to the grids according to the error function.
where the |$M_{\rm n}$| is the mass assigned to nth grid; |$z(n)$| is the position of left side of the grid, and |$z(n+1)$| is the right side; |$N(z_{\rm s},\sigma _{s})$| is the normalized Gaussian distribution.
We only consider the |$\pm 2\sigma$| range of the Gaussian profile, which provides a sufficiently large span to accurately compute the power spectrum. In cases where both the start and end points of the profile reside within the same grid, the nearest grid point (NGP) algorithm is employed to assign the source to the nearest grid regardless of the profiles.
4. Normalize the density field to a mean-centred fluctuation field; |$\delta ({\bf s})=(\rho ({\bf s})/\ \overline{\rho })-1$|
5. Use the Fast Fourier Transform (FFT) to calculate the 3D power spectrum and then spherically average the 3D power spectrum into wavebands of k to give the 1D power spectrum, and calculate the amplitude according to equation (1). If observational settings are considered, then implement the cuts of k modes in the 3D power spectrum.
6. Multiply the power spectrum of H i density fluctuations with the square of the mean brightness temperature |$\overline{T}^2_b$|.
We assign each stretched H i source in our particle-mesh code following the NGP methods. Notably, we do not adopt the cloud-in-cell (CIC) or other high-order smoothing methods (hereafter referred to as ‘CIC methods’) and only adopt NGP due to several factors. First, as long as the sampling frequency is much larger than the frequency of our concern, NGP can be accurate enough. Secondly, as some low-pass filters are narrower than NGP (e.g. see fig. 8 in Cunnington & Wolz 2024), CIC methods damp the power spectra at small scales and also distort the shot noise plateau if not properly deconvolved. Furthermore, Fourier aliasing effects become too pronounced if the CIC or other kernels are deconvolved without performing the interlacing method to remove the aliasing effects. Hence, although the NGP methods can slightly suppress the power spectra in the smallest scales, given the limited computation cost, we believe that NGP is the optimal compromise.
3.3 Emission line profile model
The emission line profiles are generally determined both by the velocity profile and the inclination angle i between the LOS direction and the normal to the galaxy disc. We consider two factors of the velocity profile: (1) the shape of the profile and (2) the profile width. We verify that the shape of the profile has a subtle impact on the power spectrum, while the width plays a key role.
3.3.1 Shape of emission line profile
There are two common profile shapes which correspond to two different kinematic scenarios. (1) Double Horn, found in spiral galaxies where H i is usually distributed in the outskirts of galaxies and is in regular rotation motions in the discs. Typically, these sources have a relatively large inclination angle, close to edge-on; (2) Gaussian-like distribution, found in the elliptical galaxies where the random motions in the discs dominate or in face-on H i discs. Both shapes are observed in extragalactic H i surveys (e.g. see fig. 4 in Masters et al. 2019). Beyond these two, there are also some other irregular shapes that cannot be well defined. To explore the impact of the shape of velocity profiles on the H i power spectrum, we use three toy models corresponding to three shapes, Gaussian, double-horn, and flat as seen in the right panel of Fig. 2. Shown as a schematic, the coloured lines in the right panel provide insight into how these shapes look like. The colours of the spectra lines in the left and middle panels correspond to the shapes depicted in matching colours in the right panel: pink for the Gaussian distribution, blue for uniform distribution, and yellow for double-horn distribution. The double-horn function is the same as the fiducial model of Sarkar & Bharadwaj (2019) (see equation 5 in this paper); the |$\sigma$| parameter in the Gaussian function is determined by the assigned |$W_{\rm 50}$|.
H i power spectrum in redshift space with different functions of velocity profiles. As illustrated in the right panel, different colours represent different functions of the profiles: pink for Gaussian distribution, blue for uniform distribution, and the yellow for Double Horn distribution. Especially, the grey lines correspond to the case without assigning any profile to the H i sources and are shown only for comparison. To be general, two redshifts |$z=0$| (left) and |$z = 1$| (mid) are displayed in the figure.
The left and middle panels of Fig. 2 illustrate how different shapes affect H i power spectrum at two example redshifts, |$z=0$| (left panel) and |$z=1$| (mid panel). For a fair comparison, each H i source is assigned a fixed profile width, denoted as |$W_{\rm 50}$|, i.e. the FWHM, determined by an empirical model (see more details of this part in Section 3.3.2). The mass of each H i source is equal to the integral of the profile, which means that the amplitudes of the profiles are proportional to the corresponding H i mass.
The differences between the coloured lines in the left and middle panels are subtle, which means that the impact of the profile shape on the H i power spectrum can be neglected. We also show the power spectrum of H i sources without the velocity profile as grey lines in all panels. The large differences between the coloured lines and the grey lines emphasize the importance of including the profile width |$W_{\rm 50}$| for small scales. For large scales with |$k \lt 0.5\ {\rm Mpc^{-1}}$|, the power spectra are not affected at these redshifts.
3.3.2 Emission line width |$W_{\rm 50}$|
|$W_{\rm 50}$|, also known as FWHM, is a common measure used to describe the width of emission line profiles in extragalactic H i observations (e.g. see Barnes et al. 2001; Zavala et al. 2009; Zwaan, Meyer & Staveley-Smith 2010; Haynes et al. 2018; Oman 2022; Zhang et al. 2024). By comparing the grey lines and the coloured lines in Fig. 2, it can be inferred that the value of |$W_{\rm 50}$| significantly influences the H i power spectrum at small scales. Since |$W_{\rm 50}$| is not available in the NUM catalogue, we have to quantitatively determine the |$W_{\rm 50}$| for each H i disc.
According to the Newtonian dynamics, the cold gas skirts surrounding spiral galaxies should share the same gravity potential with their resident collision-less cold dark matter (CDM) halo. Therefore, it is commonly assumed that the velocities measured in the flat part of the rotation curves |$v_c$| of spiral galaxies and the surrounding gas are equal to the maximum velocity |$V_{\rm max}$| of their resident (sub)haloes (e.g. see Gonzalez et al. 2000; Zavala et al. 2009; Chauhan et al. 2019). For these H i-rich galaxies, the emission line width |$W_{\rm 50}$| can be approximated as
where i is the random inclination angle between our LOS and the norm of disc, assumed to be entirely random between 0 and |$\pi$|, not perfect, but economical and stable. Although elliptical galaxies introduce additional complexity, this work ignores the difference between them and spiral galaxies due to their limited H i content.
Under this assumption, it is supposed that the velocity function, representing the abundance of galaxies as a function of their circular velocity |$v_c$|, follows the prediction from the CDM model (e.g. Cole & Kaiser 1989; Gonzalez et al. 2000; McGaugh 2012). However, discrepancies have arisen in the low-velocity end. Measurements from the ALFALFA survey suggest that there is a serious deficiency for galaxies at the low-velocity end compared to what is expected from CDM simulations (Zavala et al. 2009; Klypin et al. 2015). Several hypotheses have been proposed to address this problem. First, the warm dark matter can suppress the abundance of small dark matter haloes since it is more difficult for small haloes to capture WDM particles with larger velocity and smaller mass compared with CDM particles (Bode, Ostriker & Turok 2001; Zavala et al. 2009; Klypin et al. 2015); Secondly, the limited SNRs and the selection effects in the H i extragalactic surveys are also considered as ways out (e.g. Chauhan et al. 2019; Oman 2022; Sardone et al. 2024); Thirdly, it is also possible that for low-velocity galaxies/H i sources, the equal velocity assumption might fail, i.e. |$v_c \ne v_{\rm max}$| due to various baryonic effects (e.g. see Macciò et al. 2016; Brooks et al. 2017; Verbeke et al. 2017). Given the ongoing debates surrounding the VF function, this work adopts the |$v_c \ne v_{\rm max}$| assumption and expresses |$W_{\rm 50}$| as a double power-law function of |$v_{\rm max}$|.
where i is the random inclination angle between our LOS and the norm of disc, and |${A_{W_{\rm 50}}, V_0, \alpha , \beta }$| are the free parameters in this empirical model.
Fig. 3 displays the width function of |$W_{\rm 50}$| obtained from ALFALFA surveys (pink and blue circles with error bars) and our mock data at |$z=0$| (solid black lines). The blue circles are the measurements obtained from ALFALFA in this work and we use them to constrain the formula, while the pink circles are the measurements obtained by Oman (2022) shown here for reference. The grey dotted line and the dashed dotted line represent the frequency resolution of MeerKAT shown in Fig. 1.
Width function of |$W_{\rm 50}$| obtained from ALFALFA (circles with error bars) and mock data (black solid line). Pink circles with error bars are the results of Oman (2022), and blue circles are obtained in this work from the data from ALFALFA. The black lines are constructed by |$V_{\rm max}$| in the NUM catalogue data. For reference, grey dotted and dashed lines are the frequency resolution of MeerKAT at |$z=0.32$| and |$z=0.44$|, respectively.
We use the width function of |$W_{\rm 50}$| measured by the ALFALAF survey to constrain the parameters and find the best-fitting values as |$A_{W_{\rm 50}}=4.868$|, |$V_0=76.053$||$\mathrm{km\, s^{-1}}$|, |$\alpha =2.075$|, |$\beta =0.675$|. The ALFALFA survey focuses primarily on the local Universe (|$z\le 0.05$|), thereby lacking constraints on the width function at higher redshifts. However, we assume that |$W_{\rm 50}$| only evolves as |$V_{\rm max}$|, which is considered properly in the simulation. We employ the best-fitting parameters as our fiducial model to estimate |$W_{\rm 50}$| at all low redshifts (|$0\le z\le 1$|).
We further investigate how variations in parameters within our fiducial model influence the H i power spectrum. From left to right panels, Fig. 4 illustrates how statistics including 21 cm power spectra, median of |$W_{\rm 50}$|-|$M_{\rm HI}$| relation and width function of |$W_{\rm 50}$| change with parameters in equation (7). The black lines in all panels of Fig. 4 are our fiducial model, shown as a reference, while the blue/red lines represent parameters being halved/doubled. All four parameters are tested, but only |$A_{W_{\rm 50}}$| is presented as we find that H i power spectrum is insensitive to the change of all other parameters, i.e. |$V_{0}, \alpha , \beta$|. The specific reason for this phenomenon is explained in further detail in Section 4.4.
Impacts of parameter |$A_{W_{\rm 50}}$| in the |$W_{\rm 50}$| model on the 21 cm power spectrum. The left, middle, and right panels display the power spectra, median of the |$W_{\rm 50}$|–H i relation, and the width function of |$W_{\rm 50}$|, respectively. Blue lines represent the statistics with one parameter halved, while red lines depict those with parameters doubled.
4 RESULTS
In the analysis of this section, the cylindrical power spectrum, 1D power spectrum, theoretical shot noise term, and scale-dependent H i bias are presented without incorporating any additional observational effects. All calculations are implemented in the redshift space except for the H i bias.
4.1 2D power spectrum
In Fig. 5, we present the cylindrical power spectra, referred to as 2D power spectra (2DPS), of the H i density fluctuation field. The bottom panels as the zoom-in of the upper panels illustrate the Kaiser effect, a gravity-driven compression of the structure on large scales in redshift space, which manifests itself as elongation along the parallel direction in Fourier space.
From left to right, this figure display the 2D power spectrum of 21 cm signal emitted by H i, |$P_{\rm 21cm}(k_\parallel ,k_{\perp })$|, in redshift space at various redshifts (|$z=0.0, 1.0, 2.0$|). The bottom panels are the zoom-in of the upper panels. The FoG effect is evident in the upper panel with a large range of k-modes, which exhibits a step-like pattern, while the Kaiser effect induced by large-scale compression is very pronounced in the lower panel with a small range of k-modes.
The upper panels depict the FoG effect, prominently observed in large k modes. It is worth noting that the FoG effect is more of a plane that descends along the parallel direction, not exactly a step-like pattern shown in this contour figure due to the discrete colour bar. This pattern arises due to two primary factors. First, in the perpendicular direction, the small diameters of galactic H i sources allow the H i sources to be effectively approximated as point masses, leading to Poisson noise, also known as shot noise. Secondly, the emission line profile of H i results in a gradual decline along the parallel direction. These combined effects produce the observed step-like pattern at large k modes, which would otherwise resemble a shot noise plateau.
Our predictions for the 2D power spectrum, particularly the FoG effect on small scales (|$k\gt 1\ {\rm Mpc^{-1}}$|), closely resemble the results obtained from the hydrodynamical simulations of TNG (see top two rows of fig. 22 in Villaescusa-Navarro et al. 2018).
4.2 1D Power spectrum
Fig. 6 illustrates the evolution of the spherically averaged 1D power spectrum of H i density fluctuations (left panel) and 21 cm signals (mid panel) at different redshifts. The power spectra of H i density fluctuations and 21 cm signal are the |$P_{\rm HI}$| and |$P_{\rm 21cm}$| in equation (1), respectively.
Angular averaging 1D power spectrum of H i fluctuation in redshift space (left panel) and the corresponding 21 cm signal (mid panel). Right panel display the |$\Omega _{\rm HI}$|. Different colours represent mock powers spectrum or H i abundance at different redshifts ranging from 0.0 to 1.0. The grey error bars in the right panel represents the measurements of cosmic H i abundance.
Different colours denote redshifts ranging from 0.0 to 1.0. The power spectra of H i density fluctuation, denoted as |$P_{\rm HI}$| in the left panel, are very similar at different redshifts. However, it is worth emphasizing that all spectra are displayed on a logarithmic scale, and therefore evolution regarding H i fluctuation spectrum is present, although small (see more details in Section 4.3). A more pronounced decreasing trend is displayed in the middle panel. This trend, as also observed in the measurements in P23, can be attributed to the decline in abundance of H i over time. The amplitude of 1D power spectrum of 21 cm signal is positively correlated with the H i abundance |$\Omega _{\rm HI}$| as shown by equation (1). The decrease occurs because of physical processes such as consumption by star-forming and ionization by ultraviolet photons. As shown in the right panel of Fig. 6, the mock H i abundances |$\Omega _{\rm HI}$| evolving with time are shown as coloured points and the measurements of |$\Omega _{\rm HI}$| are shown as grey dots with error bars. The measurements are obtained from the extended H i mass function (HIMF) at local Universe obtained by the HIPASS survey (Zwaan et al. 2005), the ALFALFA survey (Martin et al. 2010) and the ALFA Ultra Deep Survey (AUDS) (Freudling et al. 2011; Hoppmann et al. 2015); stacking methods for the |$z\lesssim 0.6$| (Lah et al. 2007; Delhaize et al. 2013; Rhee et al. 2013, 2016) and the Damped Lyman|$\alpha$| (DLA) systems that arise from the background quasar absorption lines at higher redshifts (Rao, Turnshek & Nestor 2006; Neeleman et al. 2016; Rao et al. 2017). More information about H i abundance are listed in the fig. 14 of Rhee et al. (2018). The coloured mock dots are generally within the errors of the observed ones.
For comparison, we apply our method to the TNG100-1 simulation, with a box size of 75 Mpc h−1≈110 Mpc, and present the corresponding predictions for the H i power spectrum in Fig. 7. We note that the threshold mass for halo containing H i is approximately |$10^{8.5} M_{\odot }$| in NUM model and |$10^{9.5} M_{\odot }$| in TNG100-1. However, we find that H i residing in small haloes (|$M_{\rm halo}\lt 10^{9.5} M_{\odot }$|) in the NUM model has only a subtle impact on the power spectrum, so the comparison is essentially fair. We use the TNG H i catalogues and apply our fiducial |$W_{\rm 50}$| model to subhalo |$v_{\rm max}$| catalogue data in TNG to determine the profile width of each H i source. The Gaussian shape is adopted for the emission line profiles. Red lines in Fig. 7 depict the outcomes derived from the TNG galaxy catalogue, and the grey lines represent our results based on the NUM catalogue for comparison. Differences between these lines are negligible, especially on small scales, indicating that our method and predictions are robust across different models and simulations. The only caveat is that the red lines are slightly higher than the grey lines because the linear H i bias of TNG simulation is 20 per cent larger than that in NUM (see more details in Section 4.3). This phenomenon validates the idea that constraints from the shot noise is likely to break the degeneracy between |$\Omega _{\rm HI}$| and |$b_{\rm HI}$| (Chen et al. 2021).
Comparison between our results and the results obtained by TNG at |$z=0,0.5,1.0$|. The grey lines are our predictions, red lines display the results obtained by applying the same method to the catalogue data of TNG simulation, and the black lines are the power spectrum of particle H i data shown in fig. 21 in Villaescusa-Navarro et al. (2018). The differences between the red and grey lines are neglectable, and the differences between the particle results and our results are also small.
To validate the precision of our findings, we also conduct a comparison with the results obtained from the particle data from the TNG100-1 simulation (see fig. 21 in Villaescusa-Navarro et al. 2018). The black lines on the left and right panels display the power spectra of the H i particle data. Only the results of two snapshots are presented, as the power spectra provided by Villaescusa-Navarro et al. (2018) are limited to integer redshifts. The differences between the H i power spectra of the TNG particles and our results (NUM) are found to be minimal for the local Universe, |$z=0$|. This close agreement reinforces the accuracy of our model and its ability to effectively replicate the particle results.
However, at higher redshifts (|$z=1$|), discrepancies emerge, although they are not as pronounced. These discrepancies are mainly due to the larger H i bias of the TNG simulation on large scales and the presence of either sizeable H i discs or sparse but non-negligible H i gas distributed outside galaxies on small scales.
4.3 H i Bias
Fig. 8 illustrates the dependence of H i bias on scales, given by |$b_{\rm HI}(k)=\sqrt{(P_{\rm HI}-P_{\rm SN})/P_{\rm m}}$|, where |$P_{\rm SN}$| is the shot noise term; |$P_{\rm m}$| is the theoretical non-linear matter power spectrum calculated using CAMB (Lewis, Challinor & Lasenby 2000) assuming the same cosmology as the simulation (Klypin et al. 2016; Behroozi et al. 2019; Guo et al. 2023). Each colour corresponds to a different redshift, ranging from 0.0 to 1.0. The grey dotted line marks the |$k\approx 0.1\ \rm Mpc^{-1}$| threshold below which the linear H i bias is determined (Wang et al. 2021). Specific values of linear bias at different redshifts are presented in Table 1. Due to the limitation of the box size, the bias with k modes less than |$0.03\ h\,\rm Mpc^{-1}$| are excluded from the calculation and are not depicted in the figure. The average of k modes less than |$0.045\, h\,\rm Mpc^{-1}$| is calculated separately in this work, rather than using the centre of the k bins. This is because the sparsity of k modes near the scale of the simulation box can cause the true average k modes to deviate from the k bin centre according to our tests.
The dependence of H i bias on scale. Different colours represent different redshifts ranging from 0.0 to 1.0. The bias is calculated by dividing the H i power spectra by the theoretical matter power spectra and then taking the square root, as indicated by the y-axis label. The shapes of the lines are similar, with amplitudes increasing as redshifts increase. The grey dotted line denotes the k modes below which the linear bias is determined.
Values of H i Bias of our simulation at different redshifts.
| Redshift | 0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| |$b_{\rm HI}$| | 0.698 | 0.716 | 0.732 | 0.768 | 0.805 | 0.843 | 0.883 | 0.930 | 0.984 | 1.04 | 1.10 |
| Redshift | 0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| |$b_{\rm HI}$| | 0.698 | 0.716 | 0.732 | 0.768 | 0.805 | 0.843 | 0.883 | 0.930 | 0.984 | 1.04 | 1.10 |
Values of H i Bias of our simulation at different redshifts.
| Redshift | 0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| |$b_{\rm HI}$| | 0.698 | 0.716 | 0.732 | 0.768 | 0.805 | 0.843 | 0.883 | 0.930 | 0.984 | 1.04 | 1.10 |
| Redshift | 0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| |$b_{\rm HI}$| | 0.698 | 0.716 | 0.732 | 0.768 | 0.805 | 0.843 | 0.883 | 0.930 | 0.984 | 1.04 | 1.10 |
Above |$k\approx 0.1\ \rm Mpc^{-1}$|, there is a dip feature that appears at |$k\approx 1\, \rm Mpc$|, as reported by Villaescusa-Navarro et al. (2018) and Spinelli et al. (2020) as well. This dip is possibly a result of the bump feature appearing in the non-linear matter power spectrum at this scale. This would also explain the less pronounced dip feature at higher redshift (z>1) as seen in Villaescusa-Navarro et al. (2018) and Spinelli et al. (2020), since non-linear effects are less dominant at higher redshifts.
It should be noted that our H i bias estimate is consistently lower than that provided by other simulations (see e.g. Villaescusa-Navarro et al. 2018; Spinelli et al. 2020). In contrast, our H i bias for the local Universe is in line with the measurements of ALFALFA reported by Li et al. (2022), resulting in a 20 per cent lower estimate.
A more realistic mock has been generated by Guo et al. (2023) to assess the robustness of the H i clustering measurements obtained from ALFALFA (refer to the right panel of fig. 6 in their paper). Unfortunately, the mock results suggest that the observation volume of the ALFALFA survey is insufficient to yield robust projected H i clustering measurements. Fortunately, the H i catalogue of the local Universe by Five-hundred Aperture Spherical Telescope (FAST) is imminent (Zhang et al. 2024). In the near future, a more accurate estimate for H i bias can be achieved by analysing the combination of ALFALFA surveys and the FAST all sky HI survey (FASHI), providing insights into H i bias, particularly within the local Universe.
4.4 H i Shot noise
4.4.1 Theoretical shot noise and H i power spectrum without shot noise
Considering the redshift space distortion (RSD) in our model, the position of each source in the redshift space is affected by the peculiar velocity, while it is also stretched due to their inner motions, as discussed in Section 3.3.2. Theoretically, the total power spectrum in the redshift space can be separated into the power spectrum without shot noise |$P^{(s)}$| and the shot noise term |$P_{\rm SN}$|. For each part, we need to consider different RSD effects. Fig. 9 displays the total power spectrum and the divided two parts. The black lines are the total power spectrum obtained from our mock data, identical to the coloured lines in Fig. 6 at the corresponding redshifts.
Simulated H i power spectrum on the top of mock galaxy catalogue and theoretical shot noise term. The black lines are the 21 cm power spectrum at |$z=0.0, 0.5, 1.0$|. The red lines are the theoretical shot noise term. Blue lines are the simulated corrected power spectrum calculated by subtracting shot noise term (red lines) from the total power spectrum (black lines).
In the following section, we describe how the shot noise term can be calculated theoretically, as displayed in the red lines of Fig. 9. For this term, the positional information is completely unnecessary; only the H i mass distribution and the emission line profile of each source are required. Following the assumption in Section 3.3.2, we replace the Dirac delta function with Gaussian functions. Consequently, we can write the density fluctuation field as follows.
where |${\rm \delta ^{D}}$| is the Dirac function to describe the point sources and |$\delta (\vec{x})$| is the normalized fluctuation or the overdensity; |$M_{\rm t}$| is the total mass of all H i sources |$M_{\rm t}=\sum _i M_i$|; V is the volume of the simulation box. We set the z axis as the LOS direction and |${\mathcal {N}}(z_i, \sigma _i)$| as the Gaussian profile. |$\sigma _i$| for each source is determined by the simulated |$W_{\rm 50}$|, similar to equation (4), and |$\sigma _i$| is computed as
where |$l_{w50}$| is the stretching length after translating the |$W_{\rm 50}$| with units |$\mathrm{km\, s^{-1}}$| to comoving distance with units |$\rm Mpc$| or |${\rm Mpc}\,h^{-1}$|.
For simplicity, we only show the LOS direction, while for the perpendicular direction Dirac delta functions still work properly.
The shot noise is calculated as
The theoretical power spectrum without shot noise term can be calculated using the simplest model as follows (Peacock 1992; Ballinger, Peacock & Heavens 1996; Sarkar & Bharadwaj 2019; Zhang et al. 2020).
where |$\beta$| equals to |$f/b_{\rm HI}$| with f being the growth factor and |$b_{\rm HI}$| being the bias of H i; the FoG effect |$D_{\rm FoG}$| is determined both by the peculiar velocity and the elongation due to the inner motions of sources. Since it is difficult to obtain the |$D_{\rm FoG}$| theoretically, in this work, we calculated |$P^{(s)}$| by subtracting the theoretical shot noise (red lines) from the simulated total power spectrum (black lines).
In particular, the power spectra without the shot noise and the shot noise term intersect at |$k\approx 1\ \rm Mpc$|, which coincides with the scale of the third measurement point of H i auto power spectra by P23. As mentioned above, there are discrepancies between the H i bias of different simulations and also measurements. Given the k mode of the intersection, it is likely that different constant H i biases would influence the power spectra at these scales. Therefore, changes in the H i bias could impact the constraints on H i abundance and parameters in the |$W_{\rm 50}$| model. More tests using simulations with a higher H i bias are necessary to further explore this issue in the future.
4.4.2 Lower limit of HIMF and the shot noise term
Taking into account the tight relation between |$v_{\rm max}$| and |$W_{\rm 50}$|, we give an empirical formula for this relation to determine the profile width for H i sources in our mock. However, in extragalactic H i surveys such as HIPASS, ALFALFA, and FASHI, it is more direct to obtain samples with combined information including H i mass and corresponding |$W_{\rm 50}$|. Given the amplifying effect of squaring on large values, it is likely to reproduce the shot noise term even if sources with small H i mass are not included.
We test how the lower limit of HIMF can impact the shot noise term using the theoretical method as equation (11) and present the results in Fig. 10. Different colours represent different lower limits of the H i mass. According to the middle panel, the knowledge of the most massive H i sources (|$M_{\rm HI} \gt 10^{9.2}\rm {\rm M}_{\odot }$|) can reproduce 99 per cent of the original shot noise.
Impact of lower limits in H i mass function (H i MF) on the theoretical shot noise term of 21 cm power spectrum. Different colours of lines represents different lower limits, as shown by the right panel. The left panel presents the theoretical shot noise with different cuts in H i MF. Mid panel show the fraction of shot noise with different lower limit in total shot noise. Right panel presents the |$W_{\rm 50}$|-H i mass relation in our mock catalogue with the coloured lines as the lower limits.
In the middle panel of Fig. 10, the decreasing feature of the fractional line can be explained by the |$W_{\rm 50}$|-H i mass relation depicted in the right panel. Massive H i sources typically maintain large |$W_{\rm 50}$|, which causes suppression to occur at smaller k modes, and vice versa. Consequently, the shot noise with larger lower limit appears to be less at large k modes due to the absence of less massive H i sources.
5 COMPARISON WITH OBSERVATION
In this section, we perform a comparative analysis between our model and the measurements reported by P23 (see more details in Section 2.2). This analysis accounts for the observational settings including frequency resolution and the horizon limit, both of which significantly influence the observed 21 cm power spectrum and help explain certain features shown in the measurement. Subsequently, simple fitting results are applied to provide initial interpretations of the measured data.
5.1 Observational settings
5.1.1 Horizon limit
Currently, the main challenge of H i intensity mapping is extracting the signal from the foreground, which is several orders of magnitude brighter than the H i signal. Fortunately, foregrounds such as Galactic synchrotron emission and radio galaxies exhibit smooth frequency characteristics. With a limited number of bright radio sources in the detection patches, P23 adopted a foreground avoidance technique to separate the foreground from the signal, that is, by properly accounting for the horizon limit. According to Liu, Parsons & Trott (2014), the horizon limit is taken into account using
where the |$\theta _0$| is the angular extent of the MeerKAT dishes. By assuming an extreme case, P23 adopted a horizon limit as |$k_\parallel \gt 0.3k_{\rm \perp }$| to ensure the robustness of the measurements.
For comparison, this work also follows this limit to truncate the cylindrical 2D power spectrum to obtain the 1D power spectrum. The upper panels of Fig. 11 illustrate the impact of the horizon limit truncation on the measured 21 cm power spectrum. The black dots with error bars represent the measurements reported by P23. Pink lines are the original prediction of the 21 cm power spectrum with no observational settings considered, while green lines are the predicted 21 cm power spectrum with a horizon limit considered, as shown by the diagram in the right panel. The horizon limit truncates the large parts of 2DPS and results in an obvious suppression on large k-modes and a subtle enhancement on small k-modes.
Top: How the horizon limit affects the measured 21 cm power spectrum at |$z=0.32,\ 0.44$|. The black dots with error bars in the left and mid panels are the measurements, while the coloured lines are the predicted 21 cm power spectrum. Pink lines represent the original power spectrum, while green lines represent the power spectra that are truncated by the horizon limit. The black line labelled as ‘Horizon Limit’ in the right panel is the one adopted by P23. Bottom: How the frequency resolution affects the measured 21 cm power spectrum at |$z=0.32,\ 0.44$|. The black dots with error bars in the left and mid panels are the measurements, while the solid coloured lines are the predicted 21 cm power spectrum. The coloured regions in the right panel present the cylindrical power spectrum that corresponds to the predicted 1D power spectra lines of the same colour in left and mid panel. For comparison, predicted power spectra with larger |$\Omega _{\rm HI}$| are also shown as the dashed coloured lines.
If we consider an extreme case without emission line profiles, the truncation by the horizon limit would not suppress the mock H i power spectra on small scales. Thus, we suggest that the strength of this suppression depends not only on the size of the truncation, but also on how we model the |$W_{\rm 50}$|.
5.1.2 Frequency resolution
Frequency resolution, given by the channel width of the receiver, is another observational effect that influences the measurements. For 21 cm surveys, the frequency of the emission line, reveals the redshift of the source, so the frequency resolution represents the space resolution in the light-of-sight direction and therefore determines the largest |$k_{\rm \parallel }$| that can be resolved as follows
where the |$\Delta L$| is the spatial resolution corresponding to the frequency resolution |$209\, \rm kHz$|; |$\Delta \nu$| is the frequency resolution, |$\nu _{\rm HI}=1420\ \mathrm{MHz}$| is the frequency of 21 cm signal; |$h=H_0\ ({\rm km\,s^{-1}\, Mpc^{-1}})/100$| is the small Hubble constant; |$E(z)=\sqrt{\Omega _m(1+z)^3+\Omega _{\Lambda }}$|.
The bottom panels of Fig. 11 illustrate the impact of frequency resolution on the measured 21 cm power spectrum. The hard cut of frequency resolution in the cylindrical power spectrum shown in the right panel leads to a turning point in the corresponding 1D power spectrum in the left and mid panels. The red region in the right panel corresponds to the red part of predicted 1D power spectra in the left and mid panels, while the orange region corresponds to the orange 1DPS. The orange lines start from the third point from the bottom and help explain the suddenly slowing trend in the measurement. It can be inferred that if the frequency resolution gets higher, the turning points in the left and mid panel are pushed to the larger k modes until the combined colour lines converge to the green lines in Fig. 11.
The fiducial model, combined with observational settings, is depicted by the red plus orange solid lines. For comparison, the red plus orange dashed lines are also shown as the power spectra with double |$\Omega _{\rm HI}$|. The dashed lines closely match the measurements, indicating that while the amplitude of our models is much lower than that of the measurements, the trend of our prediction closely matches the measurements. It is worth noting that Padmanabhan et al. (2023) developed a theoretical framework based on the H i-halo model. After fitting their results with P23, they also found the expected amplitude of the model to be lower than the measured amplitude, especially for measurements at |$z=0.44$|.
5.2 Fittings results
Given the absence of |$W_{\rm 50}$| measurements beyond |$z=0.2$| and the consequent uncertainty in the |$W_{\rm 50}$| model, we treat both |$A_{W_{\rm 50}}$| and |$\Omega _{\rm HI}$| as free parameters in our mock. To further understand and translate the measurements into constrained parameters within the model, we constrain the two parameters using the measurement reported by P23 at |$z=0.44$| via the Monte Carlo Markov Chain (MCMC) sampling method. Assuming a |$\chi$|-square likelihood, we fit the model parameters using the emcee package (Foreman-Mackey et al. 2013). We employ eight random walkers for the two parameters and set 3000 steps for each walker to ensure convergence. We drop the first 91 steps. Gaussian prior is established for |$\Omega _{\rm HI}$| according to the measurements made by Rao et al. (2017) at |$z\approx 0.46$|, and the flat prior ranging from 0 to 10 is established for |$A_{W_{\rm 50}}$|. The grey (light-grey) contour in Fig. 12 illustrates a strong degeneracy between |$A_{W_{\rm 50}}$| and |$\Omega _{\rm HI}$|. The likelihood of samples in the grey contour is very close, therefore the best-fitting parameter is not shown here.
The parameter contour constrained by the MeerKAT autocorrelation measurement at |$z\approx 0.44$|. We illustrate the correlation between different parameters with 68 per cent (95 per cent) confidence level in grey (light-grey) contours. As a reference, the black solid line is obtained by fixing the parameter |$A_{W_{\rm 50}}$| at different values and calculating the corresponding best |$\Omega _{\rm HI}$|. The firebrick dashed line represents our fiducial model for |$W_{\rm 50}$|, which is obtained from local Universe measurements. The purple star represents the intersection of |$b_{\rm HI}\Omega _{\rm HI}\times 10^3$| embedded in the NUM model and the fiducial |$A_{W_{\rm 50}}$|. The blue line and shadow area denote the estimation and error bar, respectively, derived from MeerKAT x WiggleZ cross-correlation measurements at |$z\approx 0.43$|.
To ensure robustness, we exclude the first two measurement data points with large error bars reported by P23 from the MCMC analysis. According to our test (not shown here), the contour shape is barely affected even with the first two points being included. For comparison, the solid black line represents the |$\Omega _{\rm HI}$| values obtained by fixing |$A_{W_{\rm 50}}$| at different values and computing the corresponding best-fitting |$\Omega _{\rm HI}$|. Additionally, the red dashed line represents our fiducial model for the |$W_{\rm 50}$|, obtained from fittings to the |$W_{\rm 50}$| measurements in the local Universe (see Section 3.3.2).
The constant H i bias in our simulation is fixed to |$b_{\rm HI} = 0.82$| at |$z\approx 0.44$|, and we include it into the |$\Omega _{\rm HI}b_{\rm HI}$| parameter, shown as the label of the y-axis of the figure. Given the existing disagreement regarding H i bias and the possible impacts of |$b_{\rm HI}$| on the mock |$P_{\rm 21cm}$|, the |$\Omega _{\rm HI}b_{\rm HI}$| is presented to mitigate the potential discrepancies related to H i bias.
The constraint from fitting the measurements of P23 is shown as the grey contours. It is clear that |$\Omega _{\rm HI}b_{\rm HI}$| is strongly coupled with |$A_{W_{\rm 50}}$|. Without an accurate estimate of |$A_{W_{\rm 50}}$|, it is difficult to place tight constraints on |$\Omega _{\rm HI}b_{\rm HI}$|. The constraints of |$\Omega _{\rm HI}b_{\rm HI}$| seem to be larger than the fiducial value of the NUM model (shown as the purple star) by a factor of 2 at the fiducial |$A_{W_{\rm 50}}$|. It could be possibly caused by the uncertainties of the observational measurements and the underestimate of |$\Omega _{\rm HI}b_{\rm HI}$| in the NUM model by fitting the |$\Omega _{\rm HI}$| measurements in literature. Given that NUM is well constrained by data from existing extragalactic H i surveys, more measurements from both the IM surveys and extragalactic H i surveys are needed to further confirm this inconsistency and explore the underlying causes.
For comparison, we also include another observational measurements in the figure, i.e. the MeerKAT radio telescope x WiggleZ Dark Energy Survey (abbreviated as ‘MeerKAT x WiggleZ’ in the figure) cross-correlation measurements (Cunnington et al. 2023) at an effective scale of |$k_{\rm eff} \approx 0.13 h\rm Mpc^{-1}$|. The cross-correlation applied the H i data measured by the single-dish modes of MeerKAT, so it provides no constraint for the |$A_{W_{\rm 50}}$| parameter due to the large detected scale. Our constraint result aligns closely with the cross-correlation measurements, as evidenced by the proximity of our fiducial model to the intersection of the blue and black lines. This is because, in any case, we do not expect |$W_{\rm 50}$| to vary dramatically when the redshift changes only from 0.0 to 0.44.
With an increasing number of measurements of the H i power spectrum on non-linear scales, |$P_{\rm 21 cm}$| shows promise as a statistic for imposing additional constraints on both |$\Omega _{\rm HI}$| and |$W_{\rm 50}$| or other profile width quantities like |$W_{\rm 20},W_{\rm 10}$| in future studies.
6 CONCLUSION
In this paper, we present a simulated-based framework to predict the measured 21 cm power spectrum on non-linear scales and explore impacts of different factors on the power spectrum.
We validate the capability of the measured data to resolve galactic and subhalo structures both perpendicular to and along the line of sight (LOS), as depicted in Fig. 1. We outline the detailed methodology employed in this study, focusing on generating the H i density field from catalogue data and deriving corresponding 1D power spectra. We then introduce a novel approach for simulating the emission line profile for each H i source. Our analysis reveals that variations in the shape of the emission line profile have subtle effects on the H i power spectra, as illustrated in Fig. 2, while the profile width, represented by the full width at half maximum (|$W_{\rm 50}$|), emerges as a critical factor.
To determine |$W_{\rm 50}$| for each source, we propose an empirical model that fits well the width function of |$W_{\rm 50}$| measured by the ALFALFA survey. This model expresses |$W_{\rm 50}$| as a double power-law function of the maximum velocity (|$v_{\rm max}$|) of subhaloes. This choice stems from the expectation that the widths of gas velocity profiles should closely track the velocity profiles of their host subhaloes. The model involves four free parameters: |$A_{W_{\rm 50}}$|, |$V_{0}$|, |$\alpha$|, and |$\beta$|. Among these parameters, the H i power spectrum is particularly sensitive to variations in |$A_{W_{\rm 50}}$|. Throughout our analysis, we adopt the best-fitting values of these parameters as the fiducial model.
To provide an interpretation of the measurements conducted by P23, we consider the observational settings employed in their study, including the horizon limit and frequency resolution. In Fig. 11, we visualize how these settings impact the H i power spectrum. The horizon limit notably suppresses the power spectrum on small scales, while the frequency resolution introduces a turning point at the corresponding k mode. Due to the considerable uncertainty of |$\Omega _{\rm HI}$| and |$W_{\rm 50}$| at |$z=0.44$|, we treat |$A_{W_{\rm 50}}$| and |$\Omega _{\rm HI}$| as free parameters and constrain them using the measurements. Shown as Fig. 12, the constraint result reveals a strong degeneracy between these parameters. The line slope between |$A_{W_{\rm 50}}$| and |$\Omega _{\rm HI}$| remains unchanged, regardless of whether the first two measurement points are included. Remarkably, our constrained results align with measurements of |$b_{\rm HI}\Omega _{\rm HI}$| obtained through cross-correlation analyses.
ACKNOWLEDGEMENTS
We would like to thank the anonymous referee for their helpful comments; Zhaoting Chen, Wenxiu Yang, Aishrila Mazumder, Sourabh Paul, and Keith Grainge for helpful discussion; Yu Lei, Meng Yang for useful discussion regarding the H i emission line; Zipeng Liu, Jiajun Zhang, and Zhenyuan Wang for useful discussions regarding the signal processing; Tao Jing and Weicheng Zang for useful discussion regarding the MCMC sampling method; Marta Spinelli for useful discussions concerning the H i bias measurement; Jonghwan Rhee for kindly providing the H i abundance measurement data.
ZL and YM are supported by the National SKA Program of China (grant No. 2020SKA0110401) and NSFC (grant No. 11821303). LW is a UK Research and Innovation Future Leaders Fellow [grant MR/V026437/1]. HG is supported by the National SKA Program of China (grant No. 2020SKA0110100), the CAS Project for Young Scientists in Basic Research (No. YSBR-092), and the science research grants from the China Manned Space Project with NO. CMS-CSST-2021-A02. SC is supported by a UK Research and Innovation Future Leaders Fellowship grant [MR/V026437/1]. We acknowledge the Tsinghua Astrophysics High-Performance Computing platform at Tsinghua University and the Core Facility for Advanced Research Computing at the Shanghai Astronomical Observatory for providing computational and data storage resources that have contributed to the research results reported within this paper.
DATA AVAILABILITY
The catalogues used in this study will be shared on reasonable request to the corresponding authors. Additional derived data are available in the article.
