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

,


INTRODUCTION
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. 1997Zwaan et al. , 2005;;Rosenberg & Schneider 2002;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 Corresponding author: Nissim Kanekar nkanekar@ncra.tifr.res.inmasses (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−3 Mpc −3 dex −1 , 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, EAGLE; Schaye et al. 2015;Pillepich et al. 2019;Davé 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 red-shifts 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 characterise the Hi properties of different galaxy populations out to z ≈ 1.3 (e.g.Bera et al. 2019Bera et al. , 2023a,b;,b;Chowdhury et al. 2020Chowdhury et al. , 2021Chowdhury et al. , 2022a,b),b).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 Bband 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
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 an ≈ 2 sq.deg area, divided into seven GMRT pointings.The observations (carried out over 3 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 −1 ; Newman et al. 2013) in the DEEP2 DR4 catalogue, 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 9 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 −1 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 −1 in the rest frame of the galaxy 1 .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).are clearly different.This is because the DEEP2 Survey targetted galaxies for spectroscopy down to a magnitude limit of R AB = 24.1,resulting in a bias towards 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.
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 variancebased 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 −1 .The RMS noise on each stacked Hi 21 cm cube was obtained via Monte Carlo simulations.In each realisation 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 −1 , 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 realisation of the stacked Hi 21 cm cube.This procedure was repeated to produce 10 4 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 Figure 1.The redshift distributions of the galaxies in the three MB subsamples, shown in blue.The distributions have been normalised by the total number of subcubes in each subsample.The Hi 21 cm subcubes of each MB subsample were assigned weights such that the effective redshift distribution of each subsample is identical to the redshift distribution of the faint MB subsample (orange histograms).The number of galaxies in each subsample is indicated in each panel, with the number of Hi 21 cm subcubes shown in parentheses.
the average velocity-integrated Hi 21 cm line luminosity ( L Hi dV), in units of Jy Mpc 2 km s −1 , 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 4 × 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 −1 , +180 km s −1 ].The average Hi masses of the galaxies in the three M B subsamples are listed in Table 1.

The M
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 (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 ⟨log M Hi ⟩ 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, ⟨log M Hi ⟩ 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 log⟨M Hi ⟩ in multiple galaxy subsamples.The scaling relations directly obtained from stacking analyses are thus "mean" scaling relations and  16.6 ± 3.9 18.5 ± 4.2 10.2 ± 2.8 Table 1.The average properties of the galaxies in the three MB subsamples.For each MB subsample, the rows are (1) the MB range, (2) the number of Hi 21 cm subcubes, (3) the number of galaxies, (4) the average redshift, (5) the average MB, (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, MB values, and stellar masses are all weighted averages, with the weights chosen to ensure an identical redshift distribution for the three MB subsamples.are generally different from the median scaling relations obtained from measurements of individual Hi masses (Bera et al. 2022(Bera et al. , 2023a;;Chowdhury et al. 2022b).
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 ⟨log M Hi ⟩ = log⟨M Hi ⟩ − (ln 10/2)σ 2 . (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 α med = α mean − (ln 10/2)σ 2 .The slopes of the two relations are the same.

Systematic effects
We note that our measurement of the M Hi − M B relation at z ≈ 1 could be affected by incompleteness in the galaxy sample of the DEEP2 survey.We provide below a systematic exploration of possible biases in our measurement of the M Hi − M B relation due to such sample incompleteness.
First, the DEEP2 Galaxy Redshift survey targeted galaxies for spectroscopy down to a limiting magnitude R AB = 24.1.The survey selected galaxies from a photometric catalogue with a limiting magnitude that is 1.5 mag.fainter than the threshold for DEEP2 spectroscopy and yielded a redshift success rate that is approximately independent of magnitude and colour (Newman et al. 2013).As such, photometric incompleteness of the DEEP2 galaxy sample is unlikely to be an issue for our results.
Next, the R-band selection of the DEEP2 survey yields a rest-frame B-band completeness that is a function of both redshift and colour: for a fixed M B , the survey is biased toward galaxies with increasingly bluer rest-frame (U-B) colour at higher redshifts (Newman et al. 2013).In our current study, we have restricted our sample for Hi 21 cm stacking to galaxies with M B ≤ −20, the B-band magnitude at which the DEEP2 survey is complete to blue star-forming galaxies at z ≲ 1 (Willmer et al. 2006;Newman et al. 2013).This implies that the galaxy sample at z ≳ 1 used for the Hi 21 cm stacking could be biased towards bluer galaxies, which could affect our average Hi mass measurements.However, Chowdhury et al. (2022a) find that the DEEP2 galaxies in our sample at z > 1 have average stellar masses and average star formation rates (SFRs) consistent with the star-forming main-sequence relation derived from highly complete photometric surveys at these redshifts.This indicates that the average properties of the DEEP2 galaxies in our sample are representative of the general main-sequence galaxy population over 0.74 ≲ z ≲ 1.45.
Finally, we note that our choice of weights for galaxies in the three subsamples is based on the distribution of galaxies in the faintest subsample and is hence biased toward the galaxies at z ≲ 1 (see Figure 1) for which the DEEP2 sample is complete at M B ≤ −20 (Newman et al. 2013).Indeed, galaxies at z > 1 make up ≈ 33% of each M B subsample, after the redshift-based weighting, while galaxies at z > 1.2, for which the effects of incompleteness would be the highest, make up only ≈ 10% of each subsample.Our average Hi mass measurements in the three subsamples are thus likely to be robust to any potential effects of incompleteness in the higher-redshift bins.Overall, we conclude that our measurement of the M Hi − M B relation at z ≈ 1 is likely to be robust to the effects of the DEEP2 selection criteria.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 an HiMF at z ≈ 0 in excellent agreement with that obtained from the AL-FALFA 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 normalisation (ϕ * ) of the B-band luminosity function over z ≈ 0.2 − 1 can be described, respectively, by M B * = (−21.00±0.03)+(z− 0.5) × (−1.03 ± 0.08) and log ϕ * (10 −3 Mpc −3 mag −1 ) = (−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 Bband 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. 2 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).

RESULTS AND DISCUSSION
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 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 an Hi mass greater than M Hi (n >MHi ), 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 10 M ⊙ at z ≈ 1 is a factor of 4.3 +29.5  −1.9 (asymmetric errors, from the 68% confidence interval 3 ) higher than that in the local Universe.We find that the hypothesis that the number density of galaxies with M Hi > 10 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 10 M ⊙ at z ≈ 1, the number density at z ≈ 1 is a factor of 9.0 +27.2 −2.9 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.
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 10 M ⊙ (blue curve) and M Hi > 5 × 10 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 10 M ⊙ .For example, for an assumed scatter of 0.13 dex, half 3 We note that the upper error range is significantly larger than the lower error range, here and in the later estimates.This is because the Hi mass function depends non-linearly on the values of the parameters of the M Hi − M B relation, and the quoted uncertainties include the full formal errors on the slope of this relation, from Equation 3. Specifically, the number density of galaxies with large Hi masses is inversely related to the slope of the M Hi − M B relation, with the value rising rapidly as one approaches a slope of 0. Our current measurement of the slope of the M Hi − M B relation is consistent with 0 at ≈ 1.5σ significance, implying that the upper bound on the number density of massive galaxies at z ≈ 1 is not tightly constrained.
that in the local Universe, the number density of galaxies with M Hi > 10 10 M ⊙ at z ≈ 1 is a factor of 3.9 +14.7 −1.6 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 10 M ⊙ at z ≈ 1 is a factor of 5 +247 −3 higher than that in the local Universe.In both cases, the hypothesis that the number density of galaxies with M Hi > 10 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 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.
The situation is different for the highest-mass galaxies, with M Hi > 5 × 10 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 emphasise that the GMRT-CATz1 survey covers a large sky area of ≈ 2 sq.degrees, corresponding to a sky volume of 10 7 comoving Mpc 3 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  The figure shows that the number density of galaxies with MHi > 10 10 M⊙ at z ≈ 1 is a factor of ≈ 4 − 5 higher than that at z ≈ 0, independent of the assumed scatter in the MHi − MB relation.However, the same ratio for galaxies with MHi > 5 × 10 10 M⊙ is seen clearly to be sensitive to the assumed scatter in the relation.
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 emphasised 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, EA-GLE, 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 den-sity 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 ⊙ .

SUMMARY
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 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 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 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.NK acknowledges support from the Department of Science and Technology via a Swarnajayanti Fellowship (DST/SJF/PSA-01/2012-13).AC, NK, & JNC also acknowledge the Department of Atomic Energy for funding support, under project 12-R&D-TFR-5.02-0700.NK also acknowledges many discussions on 21cm stacking and related issues with Balpreet Kaur and Apurba Bera that have contributed to this paper.2013,2018,2022) 3. THE M Hi − M B RELATION AT Z ≈ 1 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 1 Throughout this work, we consistently use a flat Lambda-Cold Dark Matter 737 cosmology, with Ωm = 0.3, Ω Λ = 0.7, and H 0 = 70 km s −1 Mpc −1 .

Figure 2 .
Figure 2. The average Hi 21 cm emission signal from star-forming galaxies in different MB subsamples at z ≈ 1.The top panels show the stacked Hi 21 cm emission spectra, at a velocity resolution of 90 km s −1 , for galaxies in the three MB subsamples, bright (left panel), intermediate (middle panel), and faint (right panel).The dashed curve in each panel shows the ±1σ error on the stacked Hi 21 cm spectrum.The bottom panels show the average Hi 21 cm emission images of the same galaxies in the three MB subsamples, obtained by integrating the Hi 21 cm emission over the central velocity channels (= 180 km s −1 to +180 km s −1 ) of the stacked spectral cubes.The circle at the bottom left of each panel shows the 90-kpc spatial resolution of the images.The contour levels are at −3.0σ (dashed), +3.0σ, and +4.0σ significance; note that there are no −3σ features in the images.The average Hi 21 cm emission signals from all three subsamples are seen to be clearly detected in both the spectra and the images.

Figure 3 .
Figure 3.The MHi − MB relation at z ≈ 1.The red points show our measurements of the average Hi mass in the three MB subsamples.The blue curve and blue shaded region show, respectively, the best fit to the measurements and the 1σ confidence interval.The black dashed line indicates the "mean" MHi − MB relation at z ≈ 0 (Dénes et al. 2014).It is clear that our measurements of the average Hi mass in the three MB subsamples are consistent with the MHi − MB relation at z ≈ 0, indicating that the relation has not evolved between z ≈ 1 and z ≈ 0. We emphasise that the figure shows the "mean" MHi − MB relation at both z ≈ 1 and z ≈ 0 (see Section 3): we have used the measured scatter of 0.26 dex in the MHi − MB relation at z ≈ 0 (Dénes et al. 2014) and Equation 2 to obtain the mean MHi − MB relation at z ≈ 0.
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(Bera et al.  2022; Chowdhury et al. 2022b).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).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 log [M Hi /M ⊙ ] =(−0.156± 0.105) × (M B + 21) + (10.049 ± 0.068) .

4.
DETERMINING THE Hi MASS FUNCTION AT Z ≈ 1 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.

Figure 4 .Figure 5 .
Figure 4. [A] The HiMF of star-forming galaxies at z ≈ 1.The blue curve shows our measurement of the HiMF at z ≈ 1, obtained by combining the GMRT-CATz1 measurement of the MHi − MB relation with the B-band luminosity function at z ≈ 1 from the ALHAMBRA survey (López-Sanjuan et al. 2017); the shaded region shows the 68% confidence interval on the HiMF.[B]The number density of galaxies at z ≈ 1 with Hi mass greater than MHi (n>M Hi ), obtained by integrating the HiMF of Panel [A] from MHi to ∞, is shown in blue, with the 68% confidence interval indicated by the blue shaded region.The dashed curve in both panels shows the same quantities at z ≈ 0 fromJones et al. (2018).The figure shows that the number density of galaxies with MHi > 10 10 M⊙ is far larger at z ≈ 1 than at z ≈ 0.