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

ABSTRACT

Damped Lyman-α absorber (DLA), or HI 21cm absorber (H21A), is an important probe to model-independently measure the acceleration of spectroscopic velocity (v_S) via the Sandage–Loeb effect. Confined by the shortage of DLAs and background radio sources (BRSs) with adequate information, the detectable amount of DLAs is ambiguous in the bulk of previous work. After differing the acceleration of scale factor (⁠|$\ddot{a}$|⁠) from the first-order time derivative of spectroscopic velocity (⁠|$\dot{v}_\mathrm{S}$|⁠), we make a statistical investigation of the amount of potential DLAs in the most of this paper. Using kernel density estimation to depict general redshift distributions of BRSs, observed DLAs and a DLA detection rate with different limitations (1.4 GHz flux, HI column density, and spin temperature), we provide fitted multiGaussian expressions of the three components and their 1σ regions by bootstrap, with a proportional constant of H21As in detected DLAs, leading to the measurable number predictions of H21As for FAST, ASKAP, and SKA1-Mid in HI absorption blind survey. In our most optimistic condition (F_1.4 GHz > 10 mJy, N_HI > 2 × 10²⁰ cm⁻², and T_S> 500 K), the FAST, AKSAP, and SKA1-Mid would probe about 80, 500, and 600 H21As, respectively.

1 INTRODUCTION

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 modelling 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 (S–L) 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, van Elst & Heinesen 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-α forest (optical method) and damped Lyman-α absorbers (DLAs, radio method, our main concern in this paper). DLAs contain abundant dense HI gas (column density N_HI > 2 × 10²⁰ cm⁻²), which absorbs Lyman-α 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 ≈ 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_HI,X) from a DLA survey in SDSS DR5, and argued that HI distribution barely evolves due to similar shapes of f(N_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_HI distribution at 〈z〉 = 2.5, and found it matches the observed opacity-corrected distribution well. Braun (2012) used opacity-corrected high-resolution images of local extragalaxies to generate f(N_HI,X) at z ≈ 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_f), Kanekar et al. (2014) carefully studied the DLA’s spin temperature (T_S) distribution in low and high redshift, within and beyond our galaxy, and the correlation between itself and N_HI, C_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_HI, X) in a specific range of N_HI. Their DLA incidence was a bit lower but consistent with other results and uncertainties. Rao et al. (2017) used 70 MgII-pre-selected 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 DLA data, Bird, Garnett & Ho (2017) revised f(N_HI,X) and dN/dX, especially the latter at z > 4. Curran found that the reciprocal of T_S can trace the star formation density ψ* (Curran 2017b), and the fraction of the cold neutral medium (∫τdν/N_HI) has relation to ψ* (Curran 2017a). Grasha et al. (2020) evaluated f(N_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 multidimensional features, such as N_HI, T_S, and f(N_HI,X), which have been used in DLA amount predictions. Integrating both the linearly interpolated f(N_HI, z) and the simulated completeness function C(N_HI, T_S, v_FWHM, z) built from the First Large Absorption Survey in HI (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-α 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-α absorbers are often discovered in intergalactic medium and likely have less peculiar acceleration (Cooke 2020). With abundant Lyman-α 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 ≳ 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 semitheoretical velocity uncertainties in one epoch. Further observations for stricter constraints are still applied and prepared.

The redshift number density of potential DLAs was produced by Yu, Zhang & Pen (2014); Yu et al. (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 active galactic nuclei (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 data sets to renew the redshift number density of background radio sources (BRSs) and DLAs. Meanwhile, many HI 21cm absorption survey programmes are on schedule or undergoing, such as the FLASH (Allison et al. 2022) in ASKAP telescope, MeerKAT Absorption Line Survey (Gupta et al. 2016, 2021) in MeerKAT telescope, and Widefield ASKAP L-band Legacy All-sky Blind Survey (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 Section 2. With kernel density estimation (KDE) to depict data and bootstrap to give the 1σ errors, we study the redshift number density of BRSs in Section 3.1 and observed DLAs in Section 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 (SKA1M) in Section 3.3. Our discussion and conclusion are in Sections 4 and 5. All the calculation in this paper is based on a fiducial Planck18 ΛCDM model (Ω_K0 = Ω_R0 = 0, Ω_M0 = 0.315, H₀ = 67.4 km s⁻¹ Mpc⁻¹) (Planck Collaboration 2020).

2 COSMIC ACCELERATION

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

The expansion rate in ΛCDM is

where z is redshift, and Ω_R0, Ω_M0, Ω_K0, and Ω_Λ0 are today’s density parameters of radiation, matter, curvature, and dark energy, respectively. The redshift drift of Hubble-flow tracers is

where t₀ = t_obs is the time of observer, a₀ = a(t₀) is today scale factor, the dot in a represents the derivative with respect to t₀, and H₀ is Hubble constant.

The change of its spectroscopic radial velocity (v_S = cz) is

The velocity drift is an indicator of |$\dot{a}$| differences (⁠|$\dot{a}_0-\dot{a}(z)$|⁠). And according to equation (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_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 Δv_S in equation (3).

From Fig. 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_a0 ≈ 0.6) relies on q(z) and the latter (z_z0 ≈ 1.9) depends on 1 + z − E(z).

From the upper panel in Fig. 2, the universe expands in acceleration at z < z_a0 (⁠|$\ddot{a}$| > 0), and the universe expands in deceleration at z > z_a0 (⁠|$\ddot{a}$| < 0). But z_z0 cannot distinguish between the two states. For example, when z < z_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_z0 to |$\mathit{ z}$| = 0. Instead, it only proves that the observed point is escaping us visually during this era.

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

After distinguishing |$\dot{v}_{\rm S}$| and |$\ddot{a}$|⁠, in the following sections, we use existing radio sources and H21A 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.

3 EXPECTATION OF HI 21CM ABSORPTION SYSTEMS

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_HI), and spin temperature (T_S) of DLAs, could be expressed as

n_R(z, F) is the number density of BRS per square degree related with redshift z and observed 1.4 GHz flux density F. n_D(z, N_HI) is the number density of observed DLAs in every possible sightline directing towards a BRS with z and HI column density N_HI. κ(z, N_HI, T_S) is the DLA detection rate at a given N_HI and T_S level. η(z) is the proportion of H21As in DLAs with z. Moreover, ΔZ(z), used in Allison (2021), excludes the background galaxy near its foreground DLA within the radial velocity of 3000 km s⁻¹:

where c is the speed of light.

3.1 Prediction of BRSs

The common redshift distribution of radio sources in predicting potential-DLA detectable number is (de Zotti et al. 2010)

which came from a z-binned polynomial fitting of the Combined EIS-NVSS Survey of Radio Sources (CENSORS) data (Brookes et al. 2008). However, the polynomial will be negative when z exceeds 3.5. Besides, the expression contains radio sources with F_1.4 GHz ≤ 10 mJy, which are faint and time-expensive in a blind survey to acquire prominent absorption lines.

The CENSORS targeted the ESO Imaging Survey (EIS) Patch D, covering a 3 × 2 deg² field of view and 150 sources pre-selected from NRAO VLA Sky Survey (NVSS) (Condon et al. 1998). Rigby et al. (2011) presented a 135-sample subset of CENSORS at a F_1.4 GHz completeness of 7.2 mJy. The extra F_1.4 GHz can provide more information and limitation in the BRS distribution.

Besides, Marcos-Caballero et al. (2013) made a z-binned gamma fitting for CENSORS. When comparing the galaxy angular power spectrum via the Bayesian evidence test, the gamma fitting outperformed the polynomial one.

3.1.1 More data sets

The CENSORS provides many samples for radio-source redshift revolution. Nonetheless, can the n_R(z) from a 6 deg² survey represent the whole sky situation? Considering the varied physical conditions towards every sky direction in every matter cluster, it still deserves verification with more radio data sets. The Leiden–Berkeley Deep Survey (LBDS)–Hercules (Waddington et al. 2001) and Combined NVSS-FIRST Galaxies-4 (CoNFIG-4) (Gendre, Best & Wall 2010) covering the different sky areas are satisfied.

The LBDS–Hercules has 64 radio sources with F_1.4 GHz > 2 mJy in a 2 deg² sky, and we use the data collected by Rigby et al. (2011). The CoNFIG-4 contains 184 radio sources with F_1.4 GHz > 50 mJy in a 52 deg² sky. The additional two sets greatly improve our analysis.

First, our radio data sets were made 10 yr ago, so we update them following the NVSS catalogue. New sources belonging to the sky coverage of CoNFIG-4 are listed in Table 1.

Then, we plot the distributions of all BRSs with z ≥ 0.1 (beyond the local supercluster, but including a sample of z = 0.0963 from CoNFIG-4) in Fig. 3, and list the count of different 1.4 GHz flux levels in Table 2. The flux counts are quite imbalanced between the CENSORS and Hercules, indicating a possible sky coverage bias.

To keep samples rich enough, we use CENSORS and Hercules to estimate the BRS z probability density with F_1.4 GHz ≥ 10 mJy and all sets to the z density with F_1.4 GHz ≥ 50 mJy.

Although our CoNFIG-4 data set contains the original 184 and three extra radio sources, only 105 of them have z ≥ 0.0963. Therefore, we re-scale the sky coverage as 105/187 * 52 ≈ 29.2 deg², to counteract the lost z information.

3.1.2 New descriptions and one example

Polynomial fitting often has intractable zero-points. Specific fitting models such as gamma require more prior knowledge. And different-length z-bin would lose diverse information in the process.

Consequently, we introduce the KDE to implement a non-parametric fitting of z probability density without bins. Particularly when using a Gaussian kernel, the final distribution can be extracted by multiGaussian fitting. In the next paragraphs, we will show the data process of CoNFIG-4 with F_1.4 GHz ≥ 50 mJy as an example. For other data sets following the same procedure, we would give results directly.

The scale of abscissa would affect the profile’s height, so we fix the gap of z as 0.1 at the start.

We use Gaussian kernel KDE (KernelDensity in sklearn) to fit the CoNFIG-4 data (both z and F_1.4 GHz). KDE has a key parameter ‘bandwidth’ equivalently controlling the width of a single Gaussian count unit, which replaces a single rectangle count unit in a histogram. Narrower bandwidth remains more features, but leads to a more flexible model and less generalization capability.

For instance, we plot the fitted probability density function (PDF) with several bandwidths in the upper panel of Fig. 4 and show their rms error count from 500 rounds’ eight-fold random-shuffle cross-validation in its lower panel. The mean rms error decreases when the bandwidth diminishes, and the rms error distribution keeps stable. However, the 0.1 bandwidth at least needs a 5-Gaussian profile to fit. Thus, we accept a bandwidth ≥0.2 but as small as possible and always prefer a simpler multiGaussian fitting.

Next, we explore what statistical model covers CoNFIG-4 better. Noticing the denser low-redshift region, we use exponential, gamma, and Weibull distribution to fit the data.

We use Kolmogorov–Smirnov (K–S) test (kstest in scipy) to compare their fittings. If it returns a p value higher than our expected level of significance (α = 0.05), we can regard the two arrays from the same distribution. It shows the p values 0.0815 of the norm and 0.0573 of exponential (exp), while gamma and Weibull both have p values lower than 0.05. But the auto K–S test can not process the KDE curve. Therefore, in manual coding, we plot the cumulative distribution function (CDF) of real 0.1-z-bin data, exponential, KDE curve of bandwidth 0.2 and 0.3 in Fig. 5, and compare the maximum vertical distances (maxVDs) between real data CDF and tested CDF. Now, we define a constant

where α is the significant level and m and n are the numbers of real samples and tested samples, respectively. We gain the KDE profile from real samples, so C(0.05,104,104) ≈ 0.1883. If maxVD is less than the constant C, we can treat the two arrays from the same distribution.

From Fig. 5, the exp and two KDEs comparatively fit real data from the aspect of maxVD. We plot their PDFs in Fig. 6, where the exp and two KDEs are distinct. The gamma and Weibull are worse than the exp because their p values are lower than 0.05. Though the exp matches the peak point well, it behaves poorly at the adjacent redshifts, causing the overestimation in our most concerned low-redshift (z< 1) region.

Finally, multiGaussian fitting exports an analytic expression for the BRS integration in equation (8). To balance the informative feature and simple model, we plot the 0.2-bandwidth fitting in Fig. 7 and choose the 3-Gaussian result.

Because we mainly concern with the KDE performance, we will ignore other statistical model fittings hereafter.

3.1.3 Probability density of BRSs

According to the H21A spectroscopic observation (Geréb et al. 2015), where three BRSs in all 32 ones showed F_1.4 GHz ≤ 50 mJy (but ≥35 mJy). Thus, 50 mJy can be a proper lower limit for the present observation capability. And we regard 10 mJy as an optimistic capability for the next generation of radio telescopes or arrays.

We make KDE for different radio sets, accumulate with two flux levels (50 and 10 mJy), check with manual K–S test, average them with the weights of their number density per square degree sky (⁠|$\sigma _\mathrm{CO_{50}}$| = 3.5962, |$\sigma _\mathrm{CE_{50}}$| = 4.1667, |$\sigma _\mathrm{HE_{50}}$| = 5.5, |$\sigma _\mathrm{CE_{10}}$| = 18, |$\sigma _\mathrm{HE_{10}}$| = 12.5), and plot these in Fig. 8, where larger bandwidths relate to simpler models, but lose more features, especially merging the first two peaks, which heavily impact our estimation in the low-z region. Therefore, we fix the bandwidth as 0.2.

We make multiGaussian fitting for a 2-dimension (z-F_1.4 GHz) 0.2-bandwidth KDE curves in Fig. 9. The 3-Gaussian is our prioritized option for simplicity. But in the high-z region (z > 4) of the 10 mJy lower limit, it loses the last peak and causes an irreparable decrease compared with the 50 mJy lower limit, even considering the two area number densities. Such inconsistency obliges us to add an extra Gaussian component to solve it, and the final 4-Gaussian fitting parameters and their 1σ errors are given in Table 3.

The fourth Gaussian in two KDEs are constrained insufficiently: the a, μ, and σ approach their set borders of 0.15 (a+), 4(μ-), and 0.5 (σ+). If we loosen the borders, the 4-Gaussian tends to become the 3-Gaussian again. Actually, our intention is to derive an approximate KDE expression for numerical integration, and the 4-Gaussian fitting traces the KDE better than the 3-Gaussian, which is enough. Besides, we find that the former three components (eight parameters, excluding σ₃) are insensitive to the slight change of limitation (keeping μ₄ > 3.5), varying less than 5 per cent (σ₃ is 20 per cent) and 10 per cent in the fitted values and errors. Moreover, the curve tails (z > 4) are quite small, only covering less than 1 per cent of the total area under the curve. Therefore, we choose to use the 4-Gaussian.

The errors in Table 3 come from our used samples, and do not reflect the uncertainty including sampling, which need the bootstrap method to convey. The bootstrap randomly re-samples from the original data set (A₀) with replacement forming a number of new data sets (⁠|$\mathscr {A}_\mathrm{N}=\lbrace A_\mathrm{i}\rbrace _\mathrm{N}$|⁠, often N ≥ 1000) with the same data volume. One conducts statistical inference from the newly generated sets (⁠|$\mathscr {A}_\mathrm{N}$|⁠) assuming the A₀ and |$\mathscr {A}_\mathrm{N}$| from the almost same underlying population. Thus, it provides a more robust 1σ estimation from the potential sample variations, particularly for small sample data.

We re-sample the three radio data sets 10 000 times separately, make number-weighted 0.2-bandwidth KDE, extract its 3σ range in every abscissa (z) point, and plot them in Fig. 10. The bootstrap 1σ errors are more inclusive than the parameter-fluctuation errors in multiGaussian (Table 3), and used by us.

The above procedure gives the radio source z probability density with different F_1.4GHz levels, which needs to multiply with a number density per square degree sky σ to become number density. With number-weighted average from the three data sets, we obtain σ′₅₀ = (105²/29.2 + 25²/6 + 11²/2)/(105 + 25 + 11) ≈ 3.714 and σ′₁₀ is 16.966. Similar to Allison et al. (2022), we also search the σ in the whole NVSS catalogue with a max radio diameter filter (θ < 60 arcsec) to exclude the extended BRSs, which worsens the covering factor of DLA. Given that the NVSS covers 82 per cent of the entire sky (33830 deg²), we collect N_50 mJy = 125525 and N_10 mJy = 544130, finally adopting σ₅₀ = 3.71 and σ₁₀ = 16.09 instead of σ′₅₀ and σ′₁₀.

3.2 Forecasting of observed H21As

The recent research of DLA distribution at wide z region (0 ≲ z ≲ 5) was advanced by Rao et al. (2017). They provided an optical MgII pre-selection criterion of DLAs, used it to identify 70 DLAs at low-z region (z ≲ 1.65) from 369 MgII absorbers, and combined many high-z (z ≳ 1.65) SDSS-detected DLAs and a low-z modified value to give a global redshift number density:

Although their result is enough for the z evolution of DLAs, they smoothed out the local DLA distribution features by three wide z bins, which is inadequate for cosmic acceleration observation. So, we use their identified 64 DLAs (N_HI > 2 × 10²⁰ cm⁻²) from 369 MgII absorbers, and miss six DLAs from their full samples. First, Rao, Turnshek & Nestor (2006) presented 41 DLAs from 197 MgII absorbers, while we only find 38 DLAs with |$W_\mathrm{0}^{\lambda 2796}\gt 0.6\mathrm{\mathring{\rm A}}$|⁠. Then, Turnshek et al. (2015) selected 26 DLAs from 96 MgII absorbers, while we just extract 23 DLAs. Neeleman et al. (2016) provided extra 23 DLAs with z and N_HI. We plot all 87 DLAs in Fig. 11.

Similarly, we make 2-dimension (z-N_HI) KDE for the 87 DLAs with two lower N_HI levels, N_HI > 2 × 10²⁰ cm⁻² (⁠|$\lg N_\mathrm{HI}\gt 20.3$| for visual simplicity, although it should be |$\lg (N_\mathrm{HI}/\mathrm{cm^{-2}})\gt 20.3$|⁠) and N_HI > 6.31 × 10²⁰ cm⁻² (⁠|$\lg N_\mathrm{HI}\gt 20.8$|⁠), and use manual K–S test to compare the constant C and their maxVDs. Provided the simple profiles of KDE in Fig. 12, we choose a 3-Gaussian fitting with 0.1 bandwidth for accuracy, and list the fitted values and errors in Table 4. A 10 000-bootstrap gives the 1σ errors in Fig. 13.

With the 21cm absorption optical depth τ₂₁(ν) = −log {1 − ΔF(ν)/[C_fF(ν)]}, N_HI can be expressed as (Kanekar et al. 2014)

where T_S and C_f is spin temperature and covering factor of DLA, ΔF is the absorbed flux density, and F is the continuum flux density. Because of the low 21cm optical depth (τ₂₁ ≪ 1) for most DLAs, equation (13) can reduce to

N_HI is determined by T_S, C_f. and integrated per cent absorption. Constraints on N_HI can be decomposed into three parts.

However, according to the same article (Kanekar et al. 2014), no statistical relations are found in T_S–N_HI and T_S–C_f, while N_HI–∫τ₂₁(ν)dν has a 3.8σ-significance correlation. They find a 37-DLA data set with their z, N_HI, estimated T_S. and C_f in the foreground of compact radio quasars. Depending on this data set, we will make a 3-dimension KDE in the parameter space of z–N_HI–T_S, and limit N_HI and T_S as before. Because the N_HI–C_f relation was not explored in their paper, we do not limit C_f together.

Grasha et al. (2020) found the best T_S is 175/C_fK for low-z (z < 1) DLAs, and 100/C_f–500/C_fK is acceptable. Allison (2021) obtained 274 and 576 K in the low-z region (z < 1) in the detection and non-detection situations (2σ). Overall, we take T_S of 500 and 1000 K as two limits, with the former two |$\lg N_\mathrm{HI}$| limits (20.3, 20.8) in our estimation of |$\kappa (z,\lg N_\mathrm{HI},T_\mathrm{S})$|⁠.

Again, we plot the extra 37 DLAs in Fig. 14. Making 3-dimension KDE with T_S and |$\lg N_\mathrm{HI}$| limits and conducting a manual K–S test, we perform a 3-Gaussian fitting (considering the actual shape of detection rate curves, we use one subtract the 3-Gaussian as the results) with 0.3 bandwidth in Fig. 15 with their parameters listed in Table 5. The less bandwidth and more Gaussians are discarded to alleviate overfitting. The 3σ errors from 10000-bootstrap are offered in Fig. 16. The 20.3–500 situation in Table 5 has the most unsuccessful constraint, reflected by the widest 3σ region in Fig. 16. But considering the good performance between the KDE curve and its multiGaussian fitting, we still accept the results.

3.3 Anticipation of potential H21As

Given all the components, we can calculate the potential H21As number density [equation (8)] in Fig. 17 and its integrated number count in Fig. 18 from z = 0.1 to eliminate the objects in the local Universe.

FAST covers about 24000 deg² (Li et al. 2018) and the redshift of 0 to 0.352 (1050–1450 MHz). ASKAP covers 33000 deg² (Allison et al. 2022) and the redshift of 0.4 to 1.0 (711.5–999.5 MHz). As for the most remarkable instrument, the SKA1M covers about 30 000 deg² (Kloeckner et al. 2015; Weltman et al. 2020) and the redshift of 0 to 1. We list the predicted detection yield in Table 6 for each instrument and lower limit condition, and the 1σ errors are multiplicative result from the three-component bootstraps.

4 DISCUSSION

We emphasize 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 centre can be accurately measured by pulsar timing (Zakamska & Tremaine 2005), Gaia Data Release 3 (Gaia Collaboration 2021) and very long baseline interferometry (VLBI) (Xu, Wang & Zhao 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, Lambert & Gontier 2011; Titov et al. 2022). The possible redshift uncertainties were given in Bolejko, Wang & Lewis (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 fig. 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 multipeak 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 Fig. 10, proving its robustness sufficiently.

The multiGaussian fittings for two components of KDEs in Tables 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 data sets, 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 per cent 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 the 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 per cent 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 with redshift ‘owners’ in Fig. 19, so it is hard to regard them as statistically less luminous in F_1.4 GHz. 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 data sets can provide a more robust prediction than a single one. Like the radio source number density per square degree sky (σ) at the end of Section 3.1.3, only with three data sets, our estimated σ₅₀′ = 3.714 can approach the NVSS value σ₅₀ = 3.71, while any single set corresponds a greater deviation.

Rao et al. (2017) provided abundant DLAs at low redshift (z ≲ 1.65) where the ground-based optical telescope can not directly observe. A reasonable speculation is that their pre-selection is biased by luminosity and dust extinction, having the risk of omitting some MgII absorbers, DLAs without MgII absorption, and deriving biasedly higher Ω_DLA. However, they found (in Section 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_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 pre-selection. 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 data sets, 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_HI, which is included in Rao’s data, and we want to use 2-dimension (z and N_HI) KDE for more features. Moreover, we plot Curran’s and Rao’s data in Fig. 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 Fig. 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 Fig. 20, indicating the low-z shortage of Rao’s samples. But the single value point n_DLA(z = 0) can not change the KDE result unless we have real low-redshift DLAs to feed the algorithm.

Our N_HI–T_S constraints for DLAs and the H21A-to-DLA proportional constant η, 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_HI, T_S, and C_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 an N_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_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_HI I-type ones are more than the high-N_HI. Namely, our higher sensitivity condition in N_HI (to discover more DLAs) is to lower N_HI, while they need to increase it. Their complicated functions |$\mathscr {C}$| and |$\mathscr {F}$|⁠, both relate to N_HI, and both are not plotted with respect to N_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 equation (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 per cent AGNs are radio-loud and ignored AGN’s z distribution.

5 CONCLUSION

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 emphasize 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_R, the radio-source redshift number density per square degree via three data sets (with z and F_1.4GHz), (ii) n_D, the DLA redshift number density in the sightline through 87 DLAs from two low-redshift (z ≲ 1.65) data sets (with z, N_HI), and (iii) κ, the DLA detection rate at given N_HI–T_S limitations from an extra 37-DLA data set and η (the fraction of H21As in DLAs). (1) Introducing the KDE method, we replace the traditional description of the n_R and n_D, and propose a data-based κ function, with multiGaussian profiles (Tables 3, 4, and 5) of KDE curves. (2) Adopting the bootstrap method, we give the former three functions 1σ errors (Figs. 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 SKA1M covers 600. Although our results are predicted H21A amount, it is convenient to turn them into DLA amount by dividing the proportional constant η = 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.

Acknowledgement

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

DATA AVAILABILITY

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

References