Statistical distribution of HI 21cm intervening absorbers as potential cosmic acceleration probes

An accelerating expansion is happening in our Universe, whose acceleration has not been measured model-independently up to now. An unprecedentedly precise observational determination will provide more clues to the fundamental modeling of the expanding mechanism, such as the most recognized dark energy or many candidates of Equation of State.

Common cosmological probes offer the model-dependent acceleration of scale factor($\ddot{a}$), while the first order time derivative of spectroscopic velocity ($\dot{v}_{\mathrm{S}}$) of objects faithfully tracing the Hubble flow in real time are still beyond our reach. Proposed by Sandage (1962) and improved by Loeb (1998), the redshift drift, namely the Sandage-Loeb (SL) effect is a model-independent probe of $\dot{v}_{\mathrm{S}}$ and free of cosmic geometry, used to distinguish cosmological models(Codur & Marinoni, 2021; Mishra, 2022) and explore the inhomogeneity and anisotropy of the Universe(Buchert et al., 2022; Heinesen & Macpherson, 2022). Moresco et al. (2022) recorded more details of redshift drift, and Melia (2022) discussed the significant differences between a zero and non-zero redshift drift universe.

Direct measurement of redshift drift mainly depends on lyman-$\alpha$ forest (optical method) and Damped Lyman-$\alpha$ Absorbers (DLAs, radio method, our main concern in this paper). DLAs contain abundant dense HI gas (column density $N_{\mathrm{HI}}>2\times 10^{20}\mathrm{cm^{-2}}$) which absorbs Lyman-$\alpha$ photons (optical) and most 21cm radiation (radio) from radio sources in their local rest frame.

There are many studies about the physical properties of DLAs. Zwaan et al. (2005) analysed 355 HI 21cm line of nearby extragalaxies (z$\approx$0) from Westerbork Synthesis Radio Telescope, showing that their physical properties and DLAs incident rate are aligned with those of high-z DLAs (z>2) from quasar spectra. Prochaska & Wolfe (2009) researched HI column density distribution function in comoving redshift f($N_{\mathrm{HI}}$,X) from a DLA survey in SDSS DR5, and argued that HI distribution barely evolves due to similar shapes of f($N_{\mathrm{HI}}$,X) in redshift from 2.2 to 5.5. Using statistical sub-sample from BOSS(a part of SDSS-III) quasar spectra at z>2, Noterdaeme et al. (2012) probed the $N_{\mathrm{HI}}$ distribution at $\langle z\rangle$=2.5, and found it matches the observed opacity-correct distribution well. Braun (2012) used opacity-corrected high-resolution images of local extragalaxies to generate f($N_{\mathrm{HI}}$,X) at z$\approx$0. The distributions f and the number of equivalent DLAs show systematic decreases compared with their high-z values. By defining a measurable covering factor ($C_{\mathrm{f}}$), Kanekar et al. (2014) carefully studied the DLA’s spin temperature ($T_{\mathrm{S}}$) distribution in low and high redshift, within and beyond our galaxy, and the correlation between itself and $N_{\mathrm{HI}}$, $C_{\mathrm{f}}$ and metallicity [Z/H]. Using quasar spectra from the Hubble Space Telescope archive, Neeleman et al. (2016) conducted the first DLA blind survey at z<1.6 and analysed f($N_{\mathrm{HI}}$, X) in a specific range of $N_{\mathrm{HI}}$. Their DLA incidence was a bit lower but consistent with other results and uncertainties. Rao et al. (2017) used 70 MgII-preselected DLAs (with z<1.65) to describe DLA number density in 0<z<5 combing the z=0 modification and extra high-z DLAs. From SDSS-III DR12 DLAs data. Bird et al. (2017) revised f($N_{\mathrm{HI}}$,X) and dN/dX, especially the latter at z>4. Curran found that the reciprocal of $T_{\mathrm{S}}$ can trace the star formation density $\psi^{*}$ (Curran, 2017b), and the fraction of the cold neutral medium ($\int\tau d\nu$/$N_{\mathrm{HI}}$) has relation to $\psi^{*}$ (Curran, 2017a). Grasha et al. (2020) evaluated f($N_{\mathrm{HI}}$,X) per comoving redshift and per column density interval from detected absorptions and upper limits from non-detection regions. These works provided precious DLA samples and studied their multi-dimensional features, such as $N_{\mathrm{HI}}$, $T_{\mathrm{S}}$ and f($N_{\mathrm{HI}}$,X), which have been used in DLA amount predictions. Integrating both the linearly interpolated f($N_{\mathrm{HI}}$, $z$) and the simulated completeness function C($N_{\mathrm{HI}}$, $T_{\mathrm{S}}$, $v_{\mathrm{FWHM}}$, $z$) built from the FLASH early survey (Allison et al., 2020), they explored the posterior distribution of spin temperature (Allison, 2021) considering different observations (Sadler et al., 2020), and the possible detection number (Allison et al., 2022) of intervening and associated HI 21cm absorbers given different constraints.

The first approach to $\dot{v}_{\mathrm{S}}$, Lyman-$\alpha$ forests (LFs), was carefully investigated by Liske et al. (2008) for the next generation instrument ESO-ELT and was recently renewed by Dong et al. (2022). Optical Lyman-$\alpha$ absorbers are often discovered in intergalactic medium and likely have less peculiar acceleration (Cooke, 2020). With abundant Lyman-$\alpha$ absorption in line forests, the LFs approach gathers much attention like the Cosmic Accelerometer project (Eikenberry et al., 2019), the ACCELERATION programme (Cooke, 2020), ESPRESSO and NEID spectrographs (Chakrabarti et al., 2022). Interestingly, Esteves et al. (2021) concluded that measuring redshift drift and constraining cosmological parameters is a dilemma. However, confined by the earth’s ionosphere, ground-based observations only receive LFs photons from $z\gtrsim 1.65$ (Kloeckner et al., 2015), which corresponds to the decelerating expansion and jerk era.

The second approach, HI 21cm Absorbers (H21As), which are usually discovered in DLAs, was first attempted practically by Darling (2012), in which they gave the best constraint (until now) on redshift drift of three magnitudes larger than theoretical prediction. Long-term frequency stability was verified in GBT and their used H21As. Not all DLAs (H21As) are suitable for $\dot{v}_{\mathrm{S}}$ observations. A common classification (Curran et al., 2016) divides DLAs into two types: associated (A-type, which are near the radio sources with statistically shallower and wider absorption profiles) and intervening (I-type, which are remote from the sources with deeper and narrower profiles), implying that I-type DLAs suffer less local effects and locate in colder and quieter regions, with fewer inner collisions and outside radiation. Therefore the I-type is more ideal for the observation of cosmic acceleration, and we abbreviate it as DLA hereafter. Jiao et al. (2020) made an HI 21cm absorption spectral observation in PARKES, advocating the necessity of consecutive (decade) high-resolution spectral observations against the high-velocity uncertainty. Lu et al. (2022) made a high-accuracy HI 21cm spectral observation with FAST as a preliminary effort to obtain a snapshot of the S-L signal, introducing semi-theoretical velocity uncertainties in one epoch. Further observations for stricter constraints are still applied and prepared.

The redshift number density of potential DLAs were produced by Yu et al. (2014, 2017) and Jiao et al. (2020), where they focused on CHIME, Tianlai and FAST respectively. Zhang et al. (2021) estimated the number of HI absorption lines from the radio luminosity function of radio-load AGNs. And recently Allison et al. (2022) evaluated it as mentioned before. The HI 21cm absorption surveys have found few new DLAs (Dutta et al., 2017), and the radio surveys obtained massive samples extended to deeper view and fainter sources (Matthews et al., 2021). Therefore it necessitates checking as many datasets to renew the redshift number density of BRSs and DLAs. Meanwhile, many HI 21cm absorption survey programs are on schedule or undergoing, such as the First Large Absorption Survey in HI (FLASH) (Allison et al., 2022) in ASKAP telescope, MeerKAT Absorption Line Survey (MALS) (Gupta et al., 2016, 2021) in MeerKAT telescope, Widefield ASKAP L-band Legacy All-sky Blind surveY (WALLABY) (Koribalski et al., 2020) in ASKAP too, will flesh our understanding of 21cm DLAs. Any advance in these aspects would modify the final anticipation of cosmic acceleration experiments in the radio approach.

In this paper, the difference in definition between the measured spectral acceleration ($\dot{v}_{\mathrm{S}}$) and the real cosmic (scale factor) acceleration ($\ddot{a}$) is stressed again in sec 2. With KDE to depict data and bootstrap to give the 1$\sigma$ errors, we study the redshift number density of Background Radio Sources (BRSs) in sec 3.1 and observed DLAs in sec 3.2 where we also explore DLA detection rate function and the fraction of H21As to DLAs, and estimate the detectable amount of potential H21As for FAST, ASKAP and SKA1-Mid in sec 3.3. Our discussion and conclusion are in sec 4 and 5. All the calculation in this paper is based on a fiducial Planck18 $\Lambda$CDM model ($\Omega_{\mathrm{K0}}=\Omega_{\mathrm{R0}}=0$, $\Omega_{\mathrm{M0}}=0.315$, $H_{0}=67.4\mathrm{km\ s^{-1}\ Mpc^{-1}}$) (Planck Collaboration, 2020).

For many papers containing the necessary formulae and derivations such as Liske et al. (2008), we give a brief description with a standard $\Lambda$CDM model.

The expansion rate in $\Lambda$CDM is:

where $z$ is redshift, and $\Omega_{\mathrm{R0}},\Omega_{\mathrm{M0}},\Omega_{\mathrm{K0}},\Omega_{\mathrm{\Lambda 0}}$ are today’s density parameters of radiation, matter, curvature and dark energy respectively. The redshift drift of Hubble-flow tracers is:

where $t_{0}=t_{\mathrm{obs}}$ is the time of observer, $a_{0}=a(t_{0})$ is today scale factor, the dot in $a$ represents the derivative with respect to $t_{0}$ , and $H_{0}$ is Hubble constant.

The change of its spectroscopic radial velocity ($v_{\mathrm{S}}=cz$) is:

The velocity drift is an indicator of $\dot{a}$ differences ($\dot{a}_{0}-\dot{a}(z)$). And according to eq. 3, we can use measured $\dot{v}_{\rm S}$ to further infer $\Delta\dot{a}$ and $\ddot{a}$ model-independently.

When we explain the cosmic expansion with general relativity (GR), the recession velocity is (Davis & Lineweaver, 2001):

where $D_{\rm C}(z)=c\int_{0}^{z}\frac{dz^{\prime}}{H(z^{\prime})}$ is the comoving distance. The first order time derivative of $v_{\rm G}$ is (Jiao et al., 2020):

In order to conveniently plot $\dot{v}_{\rm G}$, from the expression of deceleration parameter $q(z)$, we could write down $\ddot{a}(z)$:

Only $\ddot{a}$ can depict the physical process of deceleration or acceleration of the universe expansion, and have a consistent value at any point from any distance in the same moment, for the scale factor is generalized to the whole universe.

The observed $\dot{z}$ (or $\dot{v}_{\mathrm{S}}$) involves a comparison within a pair of space-time spots (today-observer and past-source). For a fixed source, its $\dot{z}$ (or $\dot{v}_{\mathrm{S}}$) would change with the observer’s position selection. Thus it is just a relative measurement of the apparent velocity changes for a specific local observer. Although we do measure it, it does not represent the real acceleration of the Universe’s expansion. Moreover, through a secular (decade) observation, we can derive a time-averaged $\dot{a}$ drift ($\overline{\Delta\dot{a}}$) from $\Delta v_{\mathrm{S}}$ in eq. 3.

From Figure 1, the zero-points of $\ddot{a}$ and $\dot{z}$ (or low-redshift approximated $\dot{v}_{\mathrm{S}}$) are very different, where the former ($z_{\mathrm{a0}}\approx 0.6$) relies on $q(z)$ and the latter ($z_{\mathrm{z0}}\approx 1.9$) depends on $1+z-E(z)$.

From the upper panel in Figure 2, the universe expands in acceleration at $z<z_{\mathrm{a0}}$ ($\ddot{a}$>0), and the universe expands in deceleration at $z>z_{\mathrm{a0}}$ ($\ddot{a}$<0). But $z_{\mathrm{z0}}$ cannot distinguish between the two states. For example, when $z<z_{\mathrm{z0}}$, the observed $\dot{v}_{\mathrm{S}}\propto(\dot{a}_{\mathrm{0}}-\dot{a}_{\mathrm{z}})$ is always positive, but it does not mean the whole universe expanding acceleratingly from $z_{\mathrm{z0}}$ to $z=0$. Instead, it only proves that the observed point is escaping us visually during this era.

From the lowest panel in Figure 2, if our Universe is genuinely dominated by a $\Lambda$CDM model, we could safely research spectroscopic acceleration ($\dot{v}$) with good approximation at $z\lesssim 2$.

After distinguishing $\dot{v}_{\rm S}$ and $\ddot{a}$, in the following sections, we use existing radio sources and H21As samples to make a realistic number prediction, for the purpose that one can observe them to better constrain $\dot{v}_{\rm S}$ via the radio approach.

The redshift number density of potential H21As in every square degree of the sky with three lowest detection limitations, i.e. flux density ($F$) of BRSs, HI column density ($N_{\mathrm{HI}}$) and spin temperature ($T_{\mathrm{S}}$) of DLAs, could be expressed as

$n_{\mathrm{R}}(z,F)$ is the number density of BRS per square degree related with redshift $z$ and observed 1.4GHz flux density $F$. $n_{\mathrm{D}}(z,N_{\mathrm{HI}})$ is the number density of observed DLAs in every possible sightline directing toward a BRS with $z$ and HI column density ${N_{\mathrm{HI}}}$. $\kappa(z,N_{\mathrm{HI}},T_{\mathrm{S}})$ is the DLA detection rate at a given ${N_{\mathrm{HI}}}$ and $T_{\mathrm{S}}$ level. $\eta(z)$ is the proportion of H21As in DLAs with $z$. Moreover, $\Delta Z(z)$, used in Allison (2021), excludes the background galaxy near its foreground DLA within the radial velocity of 3000$\mathrm{km/s}$:

where $c$ is the speed of light.

We emphasise that the S-L signal is a direct measurement, and its theoretical derivation is independent of Einstein’s equation and Copernican principle (Yu et al., 2014), compared with other common cosmological probes like SNe Ia. But the different matter structures of every sightline in the signal propagation necessitates further studies of the inhomogeneity influence on the S-L signal.

The systematics of the S-L effect mainly contains two parts. First, the proper acceleration of the observer involves the earthly, solar, galactic and cosmological frames. The former two terms can be calculated in the highest level of cm/s (one-decade accumulation of S-L signal) (Wright & Eastman, 2014). The solar revolution around the galactic center can be accurately measured by pulsar timing (Zakamska & Tremaine, 2005), Gaia Data Release 3 (Gaia Collaboration, 2021) and VLBI (Xu et al., 2012; Titov & Krásná, 2018) in mm/s/yr level which is the same as S-L signal. And the proper motion of extragalaxies is accessible by VLBI (Titov et al., 2011, 2022). The possible redshift uncertainties were given in Bolejko et al. (2019) and reviewed in Lu et al. (2022). Second, the peculiar motion of absorption line systems would impact the S-L signal. Combining the various physical conditions of DLAs (Dutta et al., 2017), the Neutral Mass (NM) ratio (Cold NM:Unstable NM:Warm NM = 28:20:52) in our galaxy (Murray et al., 2018), and Figure 14 from Kanekar et al. (2014), there do exist some DLAs with sharp absorptions as ideal targets for probing S-L signal. Additionally, the data process introduces new uncertainty when signal extraction and flux calibration. It is not the primary challenge but needs to maintain the same procedure and reform with new techniques.

Though simple statistical models like exponential can determine the number densities, and KDE would decrease generalization ability. We still use KDE for some reasons, most importantly to keep the method consistent in the data process. Second, as a non-parametric estimation, KDE needs less prior knowledge about the realistic distribution. Third, these simple-peak models are inadequate for multi-peak distribution. Fourth, when filtering data with different lower limits, some samples are discarded but close to the limits. However, KDE will record the contributions of these border-near samples. Although KDE is a sensitive method suffering more fluctuations from many factors, it provides similar features (a big bulge containing two peaks and a small bulge/peak) in the two panels of Figure 10, proving its robustness sufficiently.

The multi-Gaussian fittings for 2 components of KDEs in Table 3 and 5 offer poor constraints, but the fittings quite approach the KDEs and such a defect does not hinder our estimation.

Most used radio samples are radio-loud AGNs, but other information such as morphological and spectral types was not recorded. One can better grip their redshift distribution and radio luminosity function by dividing samples into more sub-types if acquiring more data, but now they are beyond our present reach.

With different radio datasets, it is necessary to compare their observations. The radio sources in CENSORS and Hercules were preliminarily selected by pioneering radio mapping and then probed by subsequent optical spectra. Original CENSORS (Brookes et al., 2008) optically observed 143 in all 150 sources, providing a 71% spectroscopic completeness. As for the rest, the $K$-$z$ (and $I$-$z$ for one source) relation was used to estimate their redshifts. Original Hercules (Waddington et al., 2001) with 72 sources, offered 47 spectroscopic redshifts, found 10 upper limits of redshift from literature, and estimated the broad-band photometric redshifts for the rest. Rigby et al. (2011) rearranged the two sets with some published reassessments. They selected a 135-sample subset for CENSORS with 73% spectroscopic completeness, and a 64-sample subset for Hercules with 40 spectroscopic and 20 photometric redshifts. As for CoNFIG-4 samples (Gendre & Wall, 2009; Gendre et al., 2010), they were selected from NVSS and FIRST (White et al., 1997) counterplots with complete Fanaroff–Riley morphology, and reviewed for their optical information. No observation was performed in CoNFIG papers, hence their redshift information was all from literature. For the samples without redshifts in CoNFIG-4, they have roughly similar but lower flux distribution compared to redshift "owners" in Figure 19, so it is hard to regard them as statistically less luminous in $F_{\mathrm{1.4GHz}}$. But we find that most samples without redshifts do not have SDSS magnitudes recorded in the original data.

Optical selections inserted in observation unavoidably bias radio number density, but more datasets can provide a more robust prediction than a single one. Like the radio source number density per square degree sky ($\sigma$) at the end of section 3.1.3, only with three datasets, our estimated $\sigma_{50}\prime$=3.714 can approach the NVSS value $\sigma_{50}$=3.71, while any single set corresponds a greater deviation.

Rao et al. (2017) provided abundant DLAs at low-redshift ($z\lesssim 1.65$) where the ground-based optical telescope can not directly observe. A reasonable speculation is that their preselection is biased by luminosity and dust extinction, having the risk of omitting some MgII absorbers, DLAs without MgII absorption, and deriving biasedly higher $\Omega_{\mathrm{DLA}}$. However, they found (in sec 4 of that article) no evidence for these effects at $z$<1.65. And Kanekar’s conclusion(Kanekar et al., 2014) also disapproves of this selection effect in the aspect of $T_{\mathrm{S}}$ distributions.

We notice that Dutta et al. (2017); Dutta (2019) advanced an absorption-blind or galaxy-selected approach, selecting a galaxy visually close to a background quasar (Quasar-Galaxy Pair, QGP) without the prior of any absorption toward the quasar. Their method is not intuitively biased as MgII preselection. We hope their project finds more QGPs for the distribution of radio DLASs.

In addition to Rao’s and Neeleman’s I-type DLAs, there exist other DLA datasets, e.g. Curran et al. (2016); Curran (2021). Curran’s data provide 85 A-type and 56 I-type DLAs with $z$>0.1. Although Curran’s data cover wider redshift and more samples, they are collected from different observations. The decisive reason for discarding Curran’s data is they do not contain $N_{\mathrm{HI}}$ which is included in Rao’s data, and we want to use 2-dimension ($z$ and $N_{\mathrm{HI}}$) KDE for more features. Moreover, we plot Curran’s and Rao’s data in Figure 20, and Rao’s KDE profile of I-type DLAs almost envelopes Curran’s, exceptionally raising a bulge mainly because of the subset T15 (Turnshek et al., 2015) with $z$ from 0.4 to 1.0. Besides, it is normal that total Curran’s KDE profile covers Rao’s, because the former contains extra A-type DLAs which are not suitable for S-L signal. Because over half DLAs (28 in 41) overlap Rao’s samples, we do not plot Neeleman’s KDE.

The different tendencies in the low-$z$ region between our KDE-derived and Rao’s DLA number density in the left panel of Figure 13, apart from their wide $z$-bin, are mainly caused by Rao’s modified DLA number density values at z=0 (Zwaan et al., 2005; Braun, 2012). The modification is reasonable if we consider the left end of Curran’s I-type KDE, which is slightly higher than Rao’s in Figure 20, indicating the low-$z$ shortage of Rao’s samples. But the single value point $n_{\mathrm{DLA}}(z=0)$ can not change the KDE result unless we have real low redshift DLAs to feed the algorithm.

Our $N_{\mathrm{HI}}$-$T_{\mathrm{S}}$ constraints for DLAs and the H21A-to-DLA proportional constant $\eta$, mainly rely on the 37-DLA data from Kanekar et al. (2014). Hence our results are inevitably affected by small samples and Poisson errors, and need more "full" information ($z$, $N_{\mathrm{HI}}$, $T_{\mathrm{S}}$ and $C_{\mathrm{f}}$) DLAs to overcome.

The recent study of DLAs detection yield includes the remarkable work from Allison et al. (2022), where they defined a completeness function $\mathscr{C}$ with the 1000 mock spectra recovered from FLASH early survey of the GAMA 23 field (Allison et al., 2020), and used a $N_{\mathrm{HI}}$ frequency distribution function $\mathscr{F}$ linearly interpolated from z=0 (Zwaan et al., 2005) and z>2 (Bird et al., 2017), to predict the number of I-type and A-type DLAs in the ASKAP range. Future researchers will be enlightened to improve their $N_{\mathrm{HI}}$ frequency distribution. And their completeness function reminds us that more factors, such as signal-to-noise rate and blind-survey resolution (related to H21A’s line width distribution), can be considered in prediction. These observational limitations from specific facilities are more concrete, and entail our further efforts.

An interesting difference between our prediction and the Allison’s of ASKAP is, in our 37 DLAs, the low-$N_{\mathrm{HI}}$ I-type ones are more than the high-$N_{\mathrm{HI}}$. Namely, our higher sensitivity condition in $N_{\mathrm{HI}}$ (to discover more DLAs) is to lower $N_{\mathrm{HI}}$, while they need to increase it. Their complicated functions $\mathscr{C}$ and $\mathscr{F}$, both relate to $N_{\mathrm{HI}}$, and both are not plotted with respect to $N_{\mathrm{HI}}$. Thus it is hard to say which one ($\mathscr{C}$ or $\mathscr{F}$) dominates this difference.

Given more physical constraints to potential intervening H21As for FAST, our most optimistic expectation is quite small (80), around 1 or 2 order magnitude smaller than some previous work, such as 1500 (Zhang et al., 2021) from a local luminosity function estimation and 2600 (Jiao et al., 2020) for a decade CRAFTS observation (Zhang et al., 2019) with a simpler integration of eq. 8. However, our optimistic result approaches the magnitude of 100 from ALFALFA-survey estimation (Wu et al., 2015). Additionally, Zhang’s 1500 result may omit that only about 10% AGNs are radio-loud and ignored AGN’s $z$ distribution.

In this paper, we make a small but important distinction between the global and dynamical $\ddot{a}$ and the local and observed $\dot{v}_{\mathrm{S}}$ at first, and emphasise that only the $\ddot{a}$ can express an actual expansion state for the whole Universe at one certain time, but it has not been measured model-independently.

Subsequently, in the bulk of the paper, we separately explore (i) $n_{\mathrm{R}}$, the radio-source redshift number density per square degree via three datasets (with $z$ and $F_{\mathrm{1.4GHz}}$), (ii) $n_{\mathrm{D}}$, the DLA redshift number density in the sightline through 87 DLAs from two low-redshift ($z\lesssim 1.65$) datasets (with $z$, $N_{\mathrm{HI}}$), and (iii) $\kappa$, the DLA detection rate at given $N_{\mathrm{HI}}$-$T_{\mathrm{S}}$ limitations from an extra 37-DLA dataset and $\eta$ (the fraction of H21As in DLAs). (1) Introducing the KDE method, we replace the traditional description of the $n_{\mathrm{R}}$ and $n_{\mathrm{D}}$, and propose a data-based $\kappa$ function, with multi-Gaussian profiles (Table 3, 4 and 5) of KDE curves. (2) Adopting the bootstrap method, we give the former three functions 1$\sigma$ errors (Figure 10, 13 and 16). (3) We predict the potential blind-survey detection amounts (Table 6) for H21As with various limitations. At most optimistic conditions, FAST covers 80 H21As, ASKAP covers 500, and SKA1-Mid covers 600. Although our results are predicted H21A amount, it is convenient to turn them into DLA amount by dividing the proportional constant $\eta$=0.62.

The lack of BRS and DLA samples with "full" information is one of the crucial obstructions to our study. Nevertheless, with more potential H21As discovered in the future, it is a good chance to research H21As’ physical essences and environments, as well as to revise the BRSs from a new aspect. And all this progress can propel our knowledge of the Universe, and help us to comprehend its expansion history and the underlying drivers.

We thank the anonymous referee for the kind comments that help us greatly to improve this paper. We acknowledge support from National SKA Program of China (2022SKA0110202) and National Natural Science Foundation of China (grants No. 61802428, 11929301).

The data we used in this paper are all open on the Internet. Three radio datasets (Rigby et al., 2011; Gendre et al., 2010) and three DLA datasets (Rao et al., 2006; Turnshek et al., 2015; Rao et al., 2017) are available in VizieR. The rest DLAs from Neeleman et al. (2016); Kanekar et al. (2014); Curran et al. (2016); Curran (2021) are directly extracted from their articles.