The Hi Mass Function of Star-forming Galaxies at z ≈ 1

Neutral atomic hydrogen (Hi) is the primary fuel for star formation in galaxies and a critical component of a galaxy's baryonic cycle (e.g., Péroux & Howk 2020). Understanding the evolution of the Hi content of galaxies with cosmological time is thus critical for an understanding of galaxy evolution. A basic descriptor of the Hi content of galaxies at any epoch is the "Hi mass function" (HiMF), the number density of galaxies of a given Hi mass as a function of the Hi mass (e.g., Bothun 1985; Briggs 1990; Briggs & Rao 1993; Rao & Briggs 1993; Zwaan et al. 1997; Rosenberg & Schneider 2002; Zwaan et al. 2005; Hoppmann et al. 2015; Jones et al. 2018). Over the past two decades, unbiased wide-field Hi 21 cm emission surveys have shown that a Schechter function provides a good fit to the HiMF at z ≈ 0, with an exponential decline at masses above the "knee" of the mass function and a power-law dependence at low Hi masses (e.g., Zwaan et al. 2005; Jones et al. 2018). At present, the most accurate measurement of the HiMF in the local Universe is from the ALFALFA survey with the Arecibo Telescope (Haynes et al. 2018), yielding a knee Hi mass of log(M_*/M_⊙) = 9.94 ± 0.05, a low-mass power-law slope of α = − 1.25 ± 0.1, and a normalization of ϕ_* = (4.7 ± 0.8) × 10⁻³ Mpc⁻³ dex⁻¹, where the errors are dominated by systematic uncertainties due to cosmic variance and the absolute flux density scale (Jones et al. 2018).

Unfortunately, the weakness of the Hi 21 cm line has meant that little is known about the HiMF at cosmological distances. So far, direct estimates of the HiMF via unbiased Hi 21 cm surveys are restricted to z ≈ 0.1 (Hoppmann et al. 2015; Ponomareva et al. 2023). Measuring the HiMF at high redshifts is of critical importance to understanding galaxy evolution. Indeed, while different cosmological hydrodynamic simulations (e.g., SIMBA, IllustrisTNG, and EAGLE; Schaye et al. 2015; Davé et al. 2019; Pillepich et al. 2019) do reproduce the HiMF at z ≈ 0, the results of these simulations for the HiMF are very different at z ≈ 1 and z ≈ 2 (e.g., Davé et al. 2020). Measurements of the HiMF at high redshifts thus offer an avenue to distinguish between different models of galaxy evolution, and their inbuilt assumptions and subgrid physics.

The inherent weakness of the Hi 21 cm line can be overcome by using the stacking approach, in which the Hi 21 cm emission signals from a sample of galaxies with known spectroscopic redshifts are combined to obtain the average Hi mass of the sample (Zwaan 2000; Chengalur et al. 2001). Such Hi 21 cm stacking has been used to characterize the Hi properties of different galaxy populations out to z ≈ 1.3 (e.g., Bera et al. 2019; Chowdhury et al. 2020, 2021, 2022a, 2022b; Bera et al. 2023a, 2023b). Recently, Bera et al. (2022) demonstrated that Hi 21 cm stacking can also be used to determine the HiMF, by combining a measurement of the M_Hi –M_B relation at z ≈ 0.34 (obtained from Hi 21 cm stacking) with the B-band luminosity function (ϕ(M_B )) at the same redshift. This approach, an extension of that used to obtain the early estimates of the HiMF in the local Universe (see, e.g., Briggs 1990; Rao & Briggs 1993; Zwaan et al. 2001), yielded the first measurement of the HiMF at intermediate redshifts.

In this Letter, we use data from the Giant Metrewave Radio Telescope (GMRT) Cold-Hi AT z ≈ 1 (GMRT-CATz1) survey (Chowdhury et al. 2022c) to determine the M_Hi –M_B relation for star-forming galaxies at z ≈ 1. We combine this relation with measurements of the B-band luminosity function at z ≈ 1 to obtain the first measurement of the HiMF at z ≈ 1, at the end of the epoch of cosmic noon.

The GMRT-CATz1 survey used 510 hr of total time with the Band-4 receivers of the upgraded GMRT to observe three sky fields of the DEEP2 Galaxy Redshift Survey (Newman et al. 2013), covering the frequency range ≈550–830 MHz. This allowed us to carry out an Hi 21 cm emission survey of galaxies at z = 0.74–1.45, over a ≈2 deg² area, divided into seven GMRT pointings. The observations (carried out over three GMRT cycles), the data analysis, the galaxy sample, and the Hi 21 cm emission stacking procedure are described in detail in Chowdhury et al. (2022c). We provide below, for completeness, a brief summary of the sample of galaxies and the approach used for the Hi 21 cm stacking.

The main sample of the GMRT-CATz1 survey contains 11,419 blue star-forming galaxies at z = 0.74–1.45 with accurate spectroscopic redshifts (redshift accuracy ≲62 km s⁻¹; Newman et al. 2013) in the DEEP2 DR4 catalog, obtained after excluding (i) red galaxies, based on the color–magnitude relation between rest-frame U-B color and absolute B-magnitude M_B (Willmer et al. 2006), (ii) radio-loud active galactic nuclei based on their rest-frame 1.4 GHz luminosity (Condon et al. 2002), (iii) low-mass galaxies with M_* < 10⁹ M_⊙, and (iv) galaxies whose Hi 21 cm spectra were found to be affected by systematic non-Gaussian issues (Chowdhury et al. 2022c). We emphasize that the exact choices of the thresholds for the above selection criteria do not have a significant effect on our measurement of the average Hi mass of the full sample of galaxies (Chowdhury et al. 2022c). We further restrict the sample for the present analysis to galaxies with M_B ≤ − 20, the B-band completeness limit of the DEEP2 survey for blue galaxies at z ≈ 1 (Newman et al. 2013). This yields a final sample of 10,177 blue star-forming galaxies with M_B ≤ − 20 at z ≈ 0.74–1.45.

For each of the three observing cycles, the GMRT-CATz1 survey yielded a spectral cube covering ≈550–830 MHz for each pointing that was observed in the cycle (Chowdhury et al. 2022c). Each galaxy of the sample thus typically has 2–3 observations of its Hi 21 cm emission, obtained from the different observing cycles. For each galaxy of the sample, a subcube was extracted from its main spectral cube, covering a spatial extent of ±500 kpc around the galaxy position and ±1500 km s⁻¹ around its redshifted Hi 21 cm line frequency. All subcubes have a spatial resolution of 90 kpc and a velocity resolution of 90 km s⁻¹ in the rest frame of the galaxy. ¹ In all, we obtained 25,892 subcubes for the 10,177 galaxies in the present sample.

The spatial resolution of 90 proper kpc was chosen to ensure that the average Hi 21 cm emission signal from the full sample of galaxies in the GMRT-CATz1 survey is not resolved (Chowdhury et al. 2022c). We further note that Hi 21 cm source confusion does not significantly affect our average Hi mass measurements at the spatial resolution of 90 kpc (Chowdhury et al. 2022c).

3.1. The Average Hi Mass for Galaxy Subsamples

To determine the M_Hi –M_B relation at z ≈ 1, we first divided the sample of 10,177 blue galaxies at z = 0.74–1.45 into three M_B subsamples with M_B < − 21.3 (bright), −21.3 ≤ M_B < − 20.9 (intermediate), and −20.9 ≤ M_B ≤ − 20.0 (faint). The number of galaxies and the number of Hi 21 cm subcubes in each M_B subsample are listed in Table 1, while the redshift distributions of the three M_B subsamples are shown in Figure 1. The redshift distributions of the three subsamples are clearly different. This is because the DEEP2 Survey targeted galaxies for spectroscopy down to a magnitude limit of R_AB = 24.1, resulting in a bias toward galaxies with higher B-band luminosity at higher redshifts (Willmer et al. 2006; Newman et al. 2013). We corrected for this bias by using weights in the Hi 21 cm stacking such that the effective redshift distributions of the bright and intermediate M_B subsamples are identical to the redshift distribution of the faint M_B subsample. We separately stacked the Hi 21 cm subcubes of the galaxies in the three subsamples, using the above weights to ensure that the redshift distributions of the three subsamples are identical.

Zoom In Zoom Out Reset image size

Table 1. The Average Properties of the Galaxies in the Three M_B Subsamples

	Bright	Intermediate	Faint
MB Range	[−23.8, − 21.3]	[−21.3, − 20.9]	[−20.9, − 20.0]
Number of Hi 21 cm Subcubes	9369	6074	10,419
Number of Galaxies	3710	2361	4106
Average Redshift	0.97	0.97	0.97
Average MB	−21.73	−21.09	−20.50
Average M* (×109 M⊙)	23.3	9.9	4.8
Average Hi Mass (×109 M⊙)	16.6 ± 3.9	18.5 ± 4.2	10.2 ± 2.8

Note. For each M_B subsample, the rows are (1) the M_B range, (2) the number of Hi 21 cm subcubes, (3) the number of galaxies, (4) the average redshift, (5) the average M_B , (6) the average stellar mass, and (7) the average Hi mass, measured from the stacked Hi 21 cm emission spectra of Figure 2. Note that the average redshifts, Hi masses, M_B values, and stellar masses are all weighted averages, with the weights chosen to ensure an identical redshift distribution for the three M_B subsamples.

Download table as:  ASCIITypeset image

For each M_B subsample, the stacked Hi 21 cm cube was obtained by taking a weighted mean, across all Hi 21 cm subcubes in the sample, of the measured Hi 21 cm luminosity density in each spatial and velocity pixel of the individual Hi 21 cm subcubes (Chowdhury et al. 2022c). We note that the weights used during stacking were only to ensure that the redshift distribution of the subsamples is identical; no additional variance-based weights were used. Any residual spectral baselines were subtracted out by fitting a second-order polynomial to the spectrum of each spatial pixel of the stacked Hi 21 cm cube, excluding velocity channels covering the central ±250 km s⁻¹. The rms noise on each stacked Hi 21 cm cube was obtained via Monte Carlo simulations. In each realization of the Monte Carlo runs, the redshift of each DEEP2 galaxy in the subsample was shifted by a value in the range ±1500 km s⁻¹, drawn randomly from a uniform distribution. The velocity-shifted Hi 21 cm subcubes were then stacked, using a procedure identical to that followed for the "true" Hi 21 cm stack, to obtain a realization of the stacked Hi 21 cm cube. This procedure was repeated to produce 10⁴ stacked Hi 21 cm cubes, from which we measured the rms noise on each spatial and velocity pixel.

Figure 2 shows the stacked Hi 21 cm emission images and Hi 21 cm spectra for the three M_B subsamples. We obtain detections, with ≈3.5–4.4σ statistical significance, of the average Hi 21 cm emission signal from the galaxies in each of the three subsamples. For each M_B subsample, the average Hi mass of the constituent galaxies was obtained from the stacked Hi 21 cm cube using the following procedure: (i) the central velocity channels of the cube were integrated to produce an image of the Hi 21 cm emission signal, (ii) a spectrum was obtained from the stacked Hi 21 cm cube at the location of the peak luminosity density of the above Hi 21 cm image, (iii) a contiguous range of velocity channels with Hi 21 cm emission detected at ≥1.5σ significance were selected to obtain a measurement of the average velocity-integrated Hi 21 cm line luminosity (∫L_Hi dV), in units of Jy Mpc² km s⁻¹, and (iv) the average velocity-integrated Hi 21 cm line luminosity was converted to the average Hi mass of the sample via the relation M_Hi = 1.86 × 10⁴ × ∫L_Hi dV, in units of M_⊙. For all three M_B subsamples, the velocity interval with Hi 21 cm emission detected at ≥1.5σ significance over contiguous velocity channels was found to be [−180 km s⁻¹, +180 km s⁻¹]. The average Hi masses of the galaxies in the three M_B subsamples are listed in Table 1.

Zoom In Zoom Out Reset image size

3.2. The M_Hi –M_B relation at z ≈1

We fitted a power-law relation to our measurements of the average Hi mass of star-forming galaxies at z ≈ 1 in the three M_B subsamples, following the procedures of Appendix A, to obtain the M_Hi –M_B relation. The best-fit M_Hi –M_B relation at z ≈ 1 is

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}\;\left[{M}_{{\rm{H}}{\rm\small{I}}}/{M}_{\odot }\right]=(-0.156\pm 0.105)\\ \,\times \ ({M}_{B}+21)+(10.127\pm 0.068).\end{array}\end{eqnarray} \tag{ 1 }$$

In the local Universe, Hi scaling relations such as the M_Hi –M_B relation are obtained from measurements of the Hi mass in individual galaxies, fitting a relation to $\langle \;\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}}\rangle$ as a function of M_B (e.g., Dénes et al. 2014). The distribution of Hi masses in each such galaxy subsample is typically log-normal (e.g., Catinella et al. 2018). In such cases, $\langle \;\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}}\rangle$ is equal to the logarithm of the median of the Hi masses of the subsample. The scaling relations obtained in the local Universe from measurements of individual Hi masses are thus "median" scaling relations. However, in stacking analysis, scaling relations such as the M_Hi –M_B relation of Equation (1) are based on measurements of the average Hi mass and hence of $\mathrm{log}\,\langle {M}_{{\rm{H}}{\rm\small{I}}}\rangle$ in multiple galaxy subsamples. The scaling relations directly obtained from stacking analyses are thus "mean" scaling relations and are generally different from the median scaling relations obtained from measurements of individual Hi masses (Chowdhury et al. 2022b; Bera et al. 2022, 2023a).

The above difference between the median scaling relations obtained at z ≈ 0 from Hi 21 cm studies of individual galaxies and the mean scaling relations obtained via Hi 21 cm stacking analyses at high redshifts must be accounted for when comparing scaling relations at different redshifts. Assuming that the distribution of Hi masses is log-normal with a logarithmic scatter σ, one obtains

$$\begin{eqnarray}&&\langle \;\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}}\rangle =\mathrm{log}\,\langle {M}_{{\rm{H}}{\rm\small{I}}}\rangle -(\mathrm{ln}\;10/2)\,{\sigma }^{2}.\end{eqnarray} \tag{ 2 }$$

This implies that the intercept of the mean scaling relation (α_mean) obtained from a stacking analysis is related to that of the median scaling relation (α_med) obtained from individual detections via ${\alpha }_{\mathrm{med}}={\alpha }_{\mathrm{mean}}-(\mathrm{ln}\;10/2)\,{\sigma }^{2}$. The slopes of the two relations are the same.

Figure 3 shows our measurements of 〈M_Hi 〉 in the three M_B subsamples and the fit of Equation (1) to the measurements. For comparison, the figure also shows the mean M_Hi –M_B relation for late-type galaxies at z ≈ 0 (Dénes et al. 2014), where we have used the measured scatter of 0.26 dex in the median local relation (Dénes et al. 2014) to obtain the mean relation (Chowdhury et al. 2022b; Bera et al. 2022). The figure shows that our measurements of 〈M_Hi 〉 in the three M_B subsamples at z ≈ 1 are remarkably consistent with the M_Hi –M_B relation at z ≈ 0. While the M_Hi –M_B relation appears slightly flatter at z ≈ 1 than at z ≈ 0, the slope of −0.156 ± 0.105 at z ≈ 1 is formally consistent (within 2σ significance) with the slope of −0.34 ± 0.01 at z ≈ 0 (Dénes et al. 2014).

Zoom In Zoom Out Reset image size

Finally, a measurement of the HiMF from the M_Hi –M_B relation requires the median relation (Bera et al. 2022). We assume that the scatter of the M_Hi –M_B relation is the same at z = 1 as at z = 0, i.e., we assume that the scatter in the relation at z ≈ 1 to also be 0.26 dex; the effect of this assumption is discussed later. This allows us to obtain the median M_Hi –M_B relation at z ≈ 1 from the mean M_Hi –M_B relation of Equation (1); this yields

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}\;\left[{M}_{{\rm{H}}{\rm\small{I}}}/{M}_{\odot }\right]=(-0.156\pm 0.105)\\ \,\times \ ({M}_{B}+21)+(10.049\pm 0.068).\end{array}\end{eqnarray} \tag{ 3 }$$

3.3. Systematic Effects

While modern measurements of the HiMF in the local Universe have relied on wide-area optically unbiased Hi 21 cm emission surveys (e.g., Zwaan et al. 2005; Jones et al. 2018), the early estimates of the HiMF at z ≈ 0 were based on combining the M_Hi –M_B relation with measurements of the B-band luminosity function ϕ(M_B ) (e.g., Briggs 1990; Rao & Briggs 1993; Zwaan et al. 2001). These authors performed a simple transformation of variable in ϕ(M_B ), substituting M_B with M_Hi via the observed M_Hi –M_B relation, to obtain ϕ(M_Hi ). However, Bera et al. (2022) found that this straightforward approach, which ignores the scatter in the M_Hi –M_B relation, results in an underestimation of the HiMF at the high-mass end. Bera et al. (2022) further demonstrated that combining the M_Hi –M_B relation of Dénes et al. (2014) with the local B-band luminosity function and incorporating the measured scatter of the local M_Hi –M_B relation (Dénes et al. 2014) yields a HiMF at z ≈ 0 in excellent agreement with that obtained from the ALFALFA survey (Jones et al. 2018). In passing, we note that it may also be possible to derive the HiMF via other scaling relations such as the M_Hi –M_* relation. However, the measured scatter of the local M_Hi –M_B relation (0.26 dex; Dénes et al. 2014) is far lower than that of other Hi scaling relations at z ≈ 0 (e.g., 0.4 dex for the M_Hi –M_* relation; Catinella et al. 2018), making it a better choice for the purpose of estimating the HiMF at high redshifts.

To determine the HiMF at z ≈ 1, we need an estimate of the rest-frame B-band luminosity function of blue star-forming galaxies at this redshift. For this, we used the Schechter function fit to the B-band luminosity function, ϕ_B (M_B ), of blue galaxies, obtained from the ALHAMBRA survey (López-Sanjuan et al. 2017), to estimate the number density of galaxies at a given M_B . López-Sanjuan et al. (2017) find that the redshift evolution of the "knee" (${M}_{B}^{* }$) and the normalization (ϕ*) of the B-band luminosity function over z ≈ 0.2–1 can be described, respectively, by ${M}_{B}^{* }\,=\,(-21.00\pm 0.03)\,+\,(z-0.5)\,\times (-1.03\pm 0.08)$ and $\mathrm{log}\;\left[{\phi }^{* }({10}^{-3}{\mathrm{Mpc}}^{-3}{\mathrm{mag}}^{-1})\right]$ = ( − 2.51 ± 0.03) + (z − 0.5) × ( − 0.01 ± 0.03). Further, these authors assume that the slope of the luminosity function (α) does not evolve with redshift, to measure α = 1.29 ± 0.02. We used the above relations for ${M}_{B}^{* }$ and ϕ* to estimate the B-band luminosity function of star-forming galaxies at z=0.97, the mean redshift of our three M_B subsamples (see Table 1).

In passing, we emphasize that the above errors on the parameters of the B-band luminosity function from the ALHAMBRA survey include an estimate of systematic uncertainties such as survey incompleteness, cosmic variance, etc. We take into account the above uncertainties while estimating the uncertainty in our measurement of the HiMF at z ≈ 1. Further, the B-band luminosity function of the ALHAMBRA survey is in excellent agreement with multiple independent measurements of the B-band luminosity function at similar redshifts, including from the DEEP2 survey (López-Sanjuan et al. 2017).

Bera et al. (2022) used a Monte Carlo approach to determine the HiMF at z ≈ 0.35, combining their measurements of the M_Hi –M_B relation at z ≈ 0.35 with the rest-frame B-band luminosity function at this redshift, and assuming that the scatter in the M_Hi –M_B relation at z ≈ 0.35 is the same as that at z ≈ 0. We follow a slightly different, but equivalent, convolution-based approach, to combine the B-band luminosity function ϕ(M_B ) at z ≈ 1 from the ALHAMBRA survey (López-Sanjuan et al. 2017) with our measurement of the median M_Hi –M_B relation at z ≈ 1 (Equation (3)), incorporating the effect of the scatter in the M_Hi –M_B relation, to obtain the HiMF at z ≈ 1. The procedure is detailed in Appendix B. ² We again assume that the scatter in the M_Hi –M_B relation at z ≈ 1 is 0.26 dex, as measured in the local Universe (Dénes et al. 2014).

Our measurement of the HiMF for star-forming galaxies at z ≈ 1 is shown in Figure 4(A), with the 1σ uncertainty on the HiMF indicated by the shaded region. Overlaid, for comparison, is the HiMF at z ≈ 0 from the ALFALFA survey (Jones et al. 2018). It is clear from the figure that the number density of high-mass galaxies, with M_Hi > 10¹⁰ M_⊙ (i.e., higher than the knee of the local Hi mass function; Jones et al. 2018), is systematically higher at z ≈ 1 than at z ≈ 0. This can be seen more clearly in Figure 4(B), which plots the number density of galaxies with a Hi mass greater than M_Hi (${n}_{\gt {M}_{{\rm{H}}{\rm\small{I}}}}$), obtained by integrating the HiMF of Figure 4(A) from M_Hi to ∞. The figure shows that the number density of galaxies with M_Hi > 10¹⁰ M_⊙ at z ≈ 1 is a factor of ${4.3}_{-1.9}^{+29.5}$ (asymmetric errors, from the 68% confidence interval) ³ higher than that in the local Universe. We find that the hypothesis that the number density of galaxies with M_Hi > 10¹⁰ M_⊙ does not evolve between z ≈ 1 and z ≈ 0 is ruled out at ≈99.7% confidence (equivalent to ≈3σ statistical significance, for Gaussian statistics). Further, considering even more massive galaxies, with M_Hi > 5 × 10¹⁰ M_⊙ at z ≈ 1, the number density at z ≈ 1 is a factor of ${9.0}_{-2.9}^{+27.2}$ higher than that in the local Universe; the hypothesis that the number density of such galaxies does not evolve from z ≈ 1 to z ≈ 0 is also ruled out at ≈99.7% confidence.

Zoom In Zoom Out Reset image size

We note that the above results are based on the assumption that the scatter in the M_Hi –M_B relation at z ≈ 1 is the same as that at z = 0. Figure 5 shows the ratio of the number density of galaxies at z ≈ 1 to that at z ≈ 0, for M_Hi > 10¹⁰ M_⊙ (blue curve) and M_Hi > 5 × 10¹⁰ M_⊙ (orange curve), as a function of the assumed scatter in the M_Hi –M_B relation at z ≈ 1. The figure clearly shows that changing the assumed value of the scatter has little effect on the number density of galaxies for Hi masses >10¹⁰ M_⊙. For example, for an assumed scatter of 0.13 dex, half that in the local Universe, the number density of galaxies with M_Hi > 10¹⁰ M_⊙ at z ≈ 1 is a factor of ${3.9}_{-1.6}^{+14.7}$ higher than that in the local Universe. Conversely, for an assumed scatter of 0.52 dex, twice that in the local Universe, the number density of galaxies with M_Hi > 10¹⁰ M_⊙ at z ≈ 1 is a factor of ${5}_{-3}^{+247}$ higher than that in the local Universe. In both cases, the hypothesis that the number density of galaxies with M_Hi > 10¹⁰ M_⊙ does not evolve from z ≈ 1 to z ≈ 0 is ruled out at ≈99.7% confidence. Thus, for a wide range of values of the scatter in the M_Hi –M_B relation, the number density of galaxies with M_Hi > 10¹⁰ M_⊙ at z ≈ 1 is ≈4–5 times higher than that in the local Universe. Our result that the number density of galaxies with Hi masses higher than the knee of the local Hi mass function is significantly higher at z ≈ 1 than at z ≈ 0 thus appears to be a robust one.

Zoom In Zoom Out Reset image size

The situation is different for the highest-mass galaxies, with M_Hi > 5 × 10¹⁰ M_⊙, where the inferred number density at z ≈ 1 does depend critically on the assumed scatter (see the orange curve in Figure 5). If the scatter in the M_Hi –M_B relation is significantly lower than the local value of 0.26 dex, then the number density of the most massive galaxies could be similar to that in the local Universe. Conversely, if the scatter in the relation is higher than in the local Universe, the number density of the most massive galaxies at z ≈ 1 would be ≫10 times higher than that at z ≈ 0. This is an important issue for searches for individual detections of Hi 21 cm emission in galaxies, which would be most sensitive to such massive galaxies. Deeper Hi 21 cm emission surveys in the future should yield a direct estimate of the scatter of the M_Hi –M_B relation at z ≈ 1 by measuring the median M_Hi relation, via a median stacking of the Hi 21 cm emission of galaxies in each M_B subsample (e.g., Bera et al. 2022).

We also emphasize that the GMRT-CATz1 survey covers a large sky area of ≈2 deg², corresponding to a sky volume of 10⁷ comoving Mpc³ over the redshift interval z = 0.74–1.45. Chowdhury et al. (2022c) estimate that the effect of cosmic variance on measured galaxy properties over the entire redshift interval is significantly lower than 10%. This is far lower than the statistical uncertainty on our measurement of the M_Hi –M_B relation at z ≈ 1. We thus do not expect the main conclusions of this paper to be affected by cosmic variance.

It is interesting to note that Bera et al. (2022) find that the number density of galaxies at z ≈ 0.35 with high Hi masses is lower than that in the local Universe, the opposite trend to that obtained here. A possible explanation for this apparent contradiction is that the high SFR of galaxies at z ≈ 1 combined with the low net gas accretion rate at this redshift (Chowdhury et al. 2023) led to a decline in the Hi mass of massive galaxies from z ≈ 1 to z ≈ 0.35. However, this trend may have later reversed, over the last ≈4 Gyr: Bera et al. (2023b) find that the net gas accretion rate is comparable to the SFR from z ≈ 0.35 to z ≈ 0, which could yield an increased number density of galaxies with a high Hi mass at z ≈ 0 compared to that at z ≈ 0.35. Conversely, as emphasized by Bera et al. (2022), their estimate of the HiMF at z ≈ 0.35 is based on a measurement of the M_Hi –M_B relation over a small cosmic volume, implying that their measurement could be affected by cosmic variance. Wide-area measurements of the HiMF at intermediate redshifts (z ≈ 0.5) are critical to achieving a more comprehensive understanding of the detailed evolution of the HiMF over the past 8 Gyr.

Finally, Figure 3 of Davé et al. (2020) shows the expected HiMF at z ≈ 1 for the TNG100, EAGLE, EAGLE-Recal, and SIMBA hydrodynamical simulations. It is clear from this figure that TNG100, EAGLE, and EAGLE-Recal all yield a lower number density of high Hi-mass galaxies at z ≈ 1, contrary to the results obtained in this work. SIMBA thus appears to be the only current cosmological hydrodynamical simulation that produces a higher number density of high Hi-mass galaxies at z ≈ 1, in rough agreement with the measurement of the HiMF at z ≈ 1 presented here. However, even the SIMBA results for the number density of high Hi-mass galaxies at z ≈ 1 are a factor of a few lower than the values obtained here at the highest Hi masses, ≳10^10.7 M_⊙.

In this Letter, we have determined the M_Hi –M_B relation for blue star-forming galaxies at z ≈ 1, by using the GMRT-CATz1 survey and Hi 21 cm stacking to measure the average Hi masses of galaxies in three independent M_B subsamples. Our measurement of the M_Hi –M_B relation is consistent with the power-law relation at z ≈ 0 (Dénes et al. 2014), in both slope and amplitude. We used our M_Hi –M_B relation, with a scatter assumed to be the same as that at z ≈ 0, along with the B-band luminosity function of blue galaxies at z ≈ 1 from the ALHAMBRA survey to obtain the first estimate of the Hi mass function at z ≈ 1. We find statistically significant evidence that the number density of galaxies with high Hi masses is higher at z ≈ 1 than in the local Universe. Specifically, the number density of galaxies with M_Hi > 10¹⁰ M_⊙ (i.e., higher than the knee of the local Hi mass function) at z ≈ 1 is a factor of ≈4–5 higher than that at z ≈ 0, for a wide range of values of the assumed scatter in the M_Hi –M_B relation. We rule out, at ≳99.7% confidence, the hypothesis that the number density of galaxies with M_Hi > 10¹⁰ M_⊙ at z ≈ 1 is the same as that at z ≈ 0. This is the first clear evidence for a change in the HiMF at z ≲ 1. For even higher Hi masses, M_Hi > 5 × 10¹⁰ M_⊙, the number density at z ≈ 1 is ≈9 times higher than that at z ≈ 0, assuming that the scatter in the M_Hi –M_B relation at z ≈ 1 is the same as that at z ≈ 0. The high inferred number density of galaxies at z ≈ 1 with large Hi reservoirs implies that it may be possible for deep Hi 21 cm emission surveys with today's radio telescopes to obtain detections of Hi 21 cm emission in individual galaxies at z ≈ 1, and thus obtain even more direct constraints on the HiMF toward the end of the epoch of cosmic noon.

We thank the anonymous referee for comments and suggestions that have improved this paper. We thank the staff of the GMRT who have made these observations possible. The GMRT is run by the National Centre for Radio Astrophysics of the Tata Institute of Fundamental Research. N.K. acknowledges support from the Department of Science and Technology via a Swarnajayanti Fellowship (DST/SJF/PSA-01/2012-13). A.C., N.K., & J.N.C. also acknowledge the Department of Atomic Energy for funding support, under project 12-R&D-TFR-5.02-0700. N.K. also acknowledges many discussions on 21 cm stacking and related issues with Balpreet Kaur and Apurba Bera that have contributed to this paper.

Software: Astropy (Astropy Collaboration et al. 2013, 2018, 2022).

We assume that the average Hi mass of galaxies at z ≈ 1 for a given M_B follows a relation of the following form.

$$\begin{eqnarray}&&\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}}=a\,({M}_{B}+21)+c.\end{eqnarray} \tag{ A1 }$$

We fit the above relation to our measurements of the average M_Hi in the three M_B subsamples by taking into account the distribution of M_B values in each subsample as well as the weights assigned to each galaxy. First, for given values of the parameters (a,c), we use Equation (A1) to compute the average Hi mass in each of the three M_B subsamples, 〈M_Hi (a, c)〉ⁱ (for the ith M_B subsample). Next, we calculate the χ² for the given values of (a,c) as follows:

$$\begin{eqnarray}&&{\chi }^{2}\,(a,c)=\displaystyle \sum _{i=1}^{3}\,{\,\left(\displaystyle \frac{\langle {M}_{{\rm{H}}{\rm\small{I}}}{\rangle }^{i}-\langle {M}_{{\rm{H}}{\rm\small{I}}}\,(a,c){\rangle }^{i}}{{\sigma }_{{M}_{{\rm{H}}{\rm\small{I}}}\,}^{i}}\right)}^{2},\end{eqnarray} \tag{ A2 }$$

where 〈M_Hi 〉ⁱ and ${\sigma }_{{M}_{{\rm{H}}{\rm\small{I}}}\,}^{i}$ are, respectively, the measured average Hi mass and the rms uncertainty on the average Hi mass for galaxies in the ith M_B subsample. We minimize the χ² of Equation (A2) using a steepest-descent algorithm to obtain the best-fit parameters (a,c).

We emphasize that Equation (A1) is based on the observed form of the M_Hi –M_B relation (and other Hi scaling relations) in the local Universe (e.g., Dénes et al. 2014). Our current measurement of the average Hi mass in only three M_B subsamples does not allow for an independent verification of the form of the M_Hi –M_B relation at z ≈ 1. This should be possible with deeper Hi 21 cm observations in the future that would yield 〈M_Hi 〉 measurements in a larger number of M_B subsamples.

We estimate the Hi mass function by combining the M_Hi –M_B relation and the B-band luminosity function (Briggs 1990; Rao & Briggs 1993; Zwaan et al. 2001), appropriately incorporating the effect of the scatter of the M_Hi –M_B relation (Bera et al. 2022).

We again assume, based on observational results at z ≈ 0 (Dénes et al. 2014), that the median M_Hi –M_B relation has the following form,

$$\begin{eqnarray}&&\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}}=a\,{M}_{B}+b,\end{eqnarray} \tag{ B1 }$$

with a log-normal scatter (σ), such that the probability that a galaxy with a B-band luminosity M_B has a Hi mass M_Hi is given by

$$\begin{eqnarray}&&p\,(\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}}| {M}_{B})\ {{dM}}_{B}=\displaystyle \frac{a}{\sqrt{2\pi }\sigma }\;\exp \;\left[-\displaystyle \frac{{\left\{\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}}-(a\ {M}_{B}+b)\right\}}^{2}\;}{2{\sigma }^{2}}\right]\,{{dM}}_{B}.\end{eqnarray} \tag{ B2 }$$

In the hypothetical case of no scatter, i.e., σ = 0, one can simply use Equation (B1) and the B-band luminosity function, ϕ_B (M_B ), to obtain the Hi mass function, ${\phi }_{\sigma =0}\,(\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}})$, by performing a transformation of variables,

$$\begin{eqnarray}&&{\phi }_{\sigma =0}\,(\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}})\,d\,(\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}})={\phi }_{B}\,\left(\displaystyle \frac{\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}}-b}{a}\right)\,\displaystyle \frac{d\,(\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}})}{a}.\end{eqnarray} \tag{ B3 }$$

This estimate of the Hi mass function ${\phi }_{\sigma =0}\,(\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}})$ is similar to the early estimates of the HiMF at z ≈ 0 (Briggs 1990; Rao & Briggs 1993; Zwaan et al. 2001). However, in the case of a nonzero scatter in the M_Hi –M_B relation, one must take into account the actual probability, from Equation (B2), that a galaxy with B-band luminosity (M_B ) has a Hi mass of M_Hi to obtain the Hi mass function, via the following equation:

$$\begin{eqnarray}&&\phi \,(\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}})\,d\,(\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}})=\left[{\int }_{-\infty }^{\infty }\,{\phi }_{B}\,({M}_{B})\ p\,(\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}}| {M}_{B})\ d\,({M}_{B})\,\right]\,\displaystyle \frac{d\,(\mathrm{log}\;{M}_{{\rm{H}}{\rm\small{I}}})}{a}.\end{eqnarray} \tag{ B4 }$$

We use the measurement of ϕ_B (M_B ) from the ALHAMBRA survey at z ≈ 1 (López-Sanjuan et al. 2017), our measurement of the median M_Hi –M_B relation (Equation (3)) at z ≈ 1, and an assumed scatter of 0.26 dex for the M_Hi –M_B relation (i.e., the same as that in the local Universe; Dénes et al. 2014), to numerically evaluate the convolution of Equation (B4); this yields the estimate of the HiMF at z ≈ 1 shown in Figure 4(A).

The errors on the HiMF are estimated using Monte Carlo simulations. In each Monte Carlo run, the parameters a, b in the M_Hi –M_B relation and the parameters ${\phi }^{* },{M}_{B}^{* },\alpha$ for ϕ_B are randomly drawn from Gaussian probability distributions, with standard deviations set to the estimated 1σ errors on the individual parameters. The randomly drawn values of the parameters are then used to evaluate Equation (B4) to obtain independent estimates of the HiMF at z ≈ 1. We repeat the above procedure 10⁵ times to obtain an estimate of the error on our estimate of the HiMF at z ≈ 1.