bridging the gap between ionized hydrogen flowing in from the intergalactic medium and molecular hydrogen, which serves as the primary fuel for star formation (M. P. Haynes et al. 1984; N. M. McClure-Griffiths et al. 2023). H I masses in galaxies can change due to the processes of accretion, consumption/ conversion, and expulsion. Consequently, tracking the evolution of H I over cosmic time is vital for understanding galaxy evolution. The number density of galaxies as a function of H I mass (i.e., the neutral hydrogen mass function or HIMF) and its integral, the contribution of galaxies to the cosmic H I density (Ω H I ), are two key metrics for describing the distribution and abundance of neutral gas across cosmic time. Measurements of the HIMF and Ω H I enable comparisons with semianalytic and numerical models of galaxy evolution (G. Popping et al. 2014; H.-S. Kim et al. 2015; R. Davé et al. 2020).

Beyond cosmic averages, understanding H I properties in different settings can illuminate the dependence of galaxy evolution processes on environmental factors. H I is particularly sensitive to galaxy interactions, as it is affected by hydrodynamic pressure and has a more broadly extended distribution within galaxies compared to stars and other gas phases (e.g., J. C. Mihos 2001). The H I contents and distributions of galaxies reflect histories of interaction and the impacts of multiple evolutionary mechanisms, with H I surveys providing vital data for validating simulations and understanding the environmental variations of galaxy evolution (e.g., A. Chung et al. 2009;B. W. Holwerda et al. 2011; M. G. Jones et al. 2018; T. N. Reynolds et al. 2022).

The H I 21 cm line, resulting from a hyperfine transition in the ground state of a hydrogen atom, serves as an essential tool for investigating neutral hydrogen on Galactic and extragalactic scales. Over the past few decades, observations of the 21 cm line have revolutionized our understanding of H I distributions within the local Universe and beyond. Numerous H I emission line surveys, conducted with single-dish radio telescopes (M. A. Zwaan et al. 2003Zwaan et al. , 2005; A. M. Martin et al. 2010; M. G. Jones et al. 2018), have played a pivotal role in assessing the global H I properties of nearby galaxies. The two largest untargeted H I surveys, the H I Parkes All-Sky Survey (HIPASS; M. J. Meyer et al. 2004) and the Arecibo Legacy Fast ALFA (ALFALFA; R. Giovanelli et al. 2005) survey, have measured the HIMF and Ω H I in the nearby Universe. These surveys and the subsequent interferometric analysis of A. A. Ponomareva et al. (2023) have been instrumental in establishing our knowledge of the HIMF and its properties, including the contributions of various galaxy populations to the HIMF and its dependence on environment (C. M. Moorman et al. 2014;K. Said et al. 2019; M. G. Jones et al. 2020).

Previous H I surveys have predominantly focused on sources in the local (z < 0.1) Universe due to the faintness (i.e., low Einstein A coefficient) of the 21 cm line and the limited capabilities of the telescopes that observe it. Despite significant investments of observing time, so far only a small number of galaxies have been detected in H I beyond the local Universe (X. Fernández et al. 2013Fernández et al. , 2016;;B. Catinella & L. Cortese 2015; K. M. Hess et al. 2019; A. R. Gogate et al. 2020), with the most distant individual detection at z ≈ 0.42 (H. Xi et al. 2024). Only two surveys have managed to assess the HIMF beyond the z < 0.1 Universe; the Arecibo Ultra-Deep Survey (AUDS; L. Hoppmann et al. 2015;H. Xi et al. 2021) covers 0 < z < 0.16, while the Blind Ultra-Deep H I Environmental Survey (BUDHIES; Y. L. Jaffé et al. 2013; A. R. Gogate et al. 2020) has, for the first time, used direct H I detections to construct the HIMF and calculate Ω H I at z ~0.2 (A. Gogate 2022), in two volumes centered on galaxy clusters.

Although compiling a large sample of direct H I detections at higher redshifts for statistical investigations remains an elusive goal, several ambitious surveys are aiming to address this deficiency, using new facilities such as the Australian Square Kilometre Array Pathfinder (ASKAP; A. W. Hotan et al. 2021), the Five-hundred-meter Aperture Spherical radio Telescope (FAST; R. Nan 2008), and the MeerKAT array (J. L. Jonas 2009), and the phased array feed upgrade of the Westerbork Synthesis Radio Telescope known as the APERture Tile In Focus (APERTIF; E. A. K. Adams et al. 2022). This next generation of untargeted H I surveys is set to enhance our understanding of how the HIMF and Ω H I evolve over cosmic time. Surveys such as the MeerKAT International GigaHertz Tiered Extragalactic Exploration (MIGHTEE; M. Jarvis et al. 2016; N. Maddox et al. 2021) and the Deep Investigation of Neutral Gas Origins (DINGO; J. Rhee et al. 2023) are designed to probe H I in galaxies across a broad range of redshifts and environments. Building on the precedent of the deep but narrow COSMOS H I Legacy Extragalactic Survey (CHILES; X. Fernández et al. 2016; K. M. Hess et al. 2019) with the Very Large Array, the Looking At the Distant Universe with the MeerKAT Array (LADUMA; S. Blyth et al. 2018) survey now aims to explore H I in emission up to an unprecedented z ~1.4, utilizing both direct and stacked detections, and is designed to detect thousands of galaxies at its targeted depth.

Previous analyses have inferred notable discrepancies in the shape of the HIMF in the local Universe (e.g., between ALFALFA's spring and fall volumes; M. G. Jones et al. 2018), but it is unclear whether these variations reflect fundamental differences in galaxy evolution in different environments or are expected consequences of cosmic variance. As we extend our observational reach with the current generation of H I surveys, the catalogs curated at higher redshifts are expected to contain few detections compared to surveys of the local Universe, further increasing the difficulties of drawing robust conclusions and motivating the development of improved analytical methods. To fully leverage the potential of these upcoming samples and to address some of the observed discrepancies at lower redshifts, it is imperative to refine and develop new methods for deriving the HIMF.

The present paper addresses this critical need by proposing a novel approach to improve the analysis and interpretation of HIMF estimates across different redshifts. We demonstrate its application by presenting the first determination of the HIMF using observations from LADUMA. Using LADUMA's Data Release 1 (DR1) data set, we derive the HIMF over roughly the last billion years (0 < z < 0.088) and calculate an associated Ω H I for this redshift range. The paper is organized as follows.

In Section 2, we describe the LADUMA survey, our processing of the DR1 data, and our definition of a high-purity H I sample. In Section 3, we describe how we measure the HIMF, which is parameterized in terms of a P. Schechter (1976) function; the results of our analysis are presented in Section 4. We discuss the novel aspects of our method and implications for its application at higher redshift in Section 5. A summary and conclusions are presented in Section 6. Throughout this paper, we assume a flat cold dark matter (ΛCDM) cosmology with H 0 = 70 km s -1 Mpc -1 , Ω m = 0.3, and Ω Λ = 0.7. All reported values from the literature have been rescaled as necessary for consistency with this cosmology.

LADUMA is a deep untargeted 21 cm survey of a single pointing using the MeerKAT array, with an area that expands from ~0.8 deg 2 at z H I = 0 to ~5 deg 2 at z H I = 1.4 (S. Blyth et al. 2018). The LADUMA pointing encompasses the Chandra Deep Field South and is centered at 03:32:30 -28:07:57 (J2000). The LADUMA data were processed on the ilifu facility operated by the Inter-University Institute for Data Intensive Astronomy (IDIA), using the processMeerKAT pipelinefoot_1 (J. D. Collier et al. 2021) combined with custom scripts. The data used in this paper were collected in 19 nighttime tracks with MeerKAT's L-band (0.88-1.67 GHz) receivers, whose average duration was 9 hr. Data were obtained at the native 32k channel resolution (26.123 kHz) of the MeerKAT correlator and processed at 8k resolution in three independent spectral windows (SPWs) of 880-933, 960-1161, and 1304-1420 MHz, in order to avoid regions of strong radio frequency interference. This paper focuses on the data in the highest-frequency (low-z) SPW, which corresponds to 0 < z H I < 0.088.

The full width at half maximum (FWHM) of the MeerKAT primary beamfoot_2 expands from 60¢.5 at 1420 MHz to 65¢ at 1304 MHz. The outstanding sensitivity of MeerKAT allows for the detection of H I emitters well past the FWHM of the primary beam. As a result, in this work we analyze a volume extending to the full width at quarter maximum (FWQM) of the primary beam, which increases from 83¢ to 90¢ over the frequency range of the SPW. A detailed account of the data reduction will be presented in a forthcoming paper (A. Kazemi-Moridani et al. 2025, in preparation), but we summarize key steps here.

Each track is individually calibrated and imaged, with robust weighting adjusted per track to minimize sidelobes. We combine the tracks after subtracting a continuum sky model from each in the uv plane. At its central frequency, the combined data cube has synthesized beam dimensions of 8.  0 × 7.  5 with P.A. = -34 o and an rms noise of ~33 μJy beam -1 per 104.52 kHz channel. Model subtraction leaves low-level continuum residuals, especially in the vicinity of bright continuum sources, that are not visible in a singletrack image. However, as the noise integrates down with the combination of multiple tracks, the remaining continuum residuals emerge above the lower noise. We have developed a splinefoot_3 -fitting algorithm to remove these residuals in our deepest data cube, which we will refer to as pixel-based continuum subtraction from this point on. This latter stage of continuum subtraction affects our measured H I fluxes, as it models the underlying continuum emission at locations where the line sources reside. Overestimation (underestimation) of the continuum level can result in underestimation (overestimation) of the line flux, as can also occur for alternative approaches to continuum subtraction (e.g., M. J. Meyer et al. 2004;M. P. Haynes et al. 2011). For the LADUMA DR1 cubes, this effect is minimal for sources with small velocity widths, as the spline-fitting algorithm only needs to model the continuum level over a few channels. For sources with larger velocity widths, the specific choices made in implementing the pixelbased continuum subtraction algorithm have the potential to systematically affect the total measured line fluxes, albeit at a modest (typically <5%-10%) level.

Given the need for an automated source-finder for injection/ recovery tests (see Section 3.1.1 below), we make use of the Source Finding Application (SoFiA; T. Westmeier et al. 2021) package, which offers excellent efficiency, flexibility, and reliability. Our source-finding approach is refined by maximizing the fraction of synthetic sources recovered with high fidelity in the entirety of the synthetic source population. To maintain consistency between finding real sources and finding synthetic sources, we apply the same approach to the real data and use a high-purity subsample of all of the line sources detected in the low-z SPW. We note that this high-purity subsample is considerably smaller than the full LADUMA DR1 source catalog. That catalog, which will be described in detail elsewhere (A. Kazemi-Moridani et al. 2025, in preparation), includes both a core sample (for which SoFiA parameters were iteratively adjusted to match the results of unguided sourcefinding with visual and matched-filtering methods) and a supplemental sample (containing the results of visual sourcefinding guided by prior knowledge of optical redshifts, as well as a more extensive exploration of SoFiA parameter space). Sources in the supplemental sample generally have lower signal-to-noise ratios (S/N) and therefore do not figure in the high-purity subsample used in this paper. The SoFiA parameters we use in this analysis match those used to develop the core DR1 sample and are listed in Table 1 (if different from package defaults). For these parameters, in the low-z SPW data cube, SoFiA finds ~190 candidate sources, of which we judge ~140 to be real detections based on visual inspection and crossmatching with optical catalogs. On the basis of previous predictive work by H. Roberts et al. (2021), and given the low redshifts in question, we do not expect any OH megamasers (OHMs) to be present in this list.

Out of our parent sample of ~140 sources, 89 are detected with high enough reliability (>99%; see Section 3.1.1 below) to be included in the high-purity sample. By limiting our analysis to the FWQM of the primary beam, we reduce the number of included detections from 89 to 84, ensuring that only high-S/N massive sources remain. Furthermore, in order to avoid biasing the HIMF by relying on the very small volume associated with the two lowest-mass sources (i.e., sources in the lowest-mass bin; see Section 5.1 for further details), we exclude them from our analysis as well, leaving 82 sources in the final sample used to measure the HIMF. Figure 1 shows the distribution of these 82 sources (see Section 3.1.1) as a function of redshift and on the plane of the sky. Given the rarity of rich clusters (e.g., G. O. Abell 1958) and the deficiency of H I emission in the densest environments (e.g., R. Giovanelli & M. P. Haynes 1985), we expect that our H I-selected sources are primarily located in field or intermediate-density environments. Considering LADUMA's angular resolution, we also expect that spectral line confusion will not significantly affect the recovered source population (e.g., M. G. Jones et al. 2015).

The HIMF is defined as the number density of galaxies as a function of H I mass per unit comoving volume, which is usually represented as

where dN gal is the number of galaxies with H I masses falling in a logarithmic mass bin centered on M H I and lying in the comoving volume dV. Determining the intrinsic HIMF from observed number counts is a complex problem, especially when a sample is not volume-limited. To measure the HIMF, it is essential to correct for the survey sensitivity, which requires estimating the completeness of the observed sample. Completeness is typically calculated for each H I mass bin as an estimate of the fraction of all galaxies within the associated mass range that have actually been detected. The main factors affecting completeness in our data are distance, nonuniform sensitivity across our field of view due to primary beam attenuation, random orientations of individual H I emitters on the sky (affecting the observed line widths of those sources), uncertainties in flux measurements, continuum subtraction effects, and source-finding accuracy. The combined effects of these factors bias our observed H I galaxy sample toward nearby, gas-rich, and low-inclination sources near the center of the LADUMA field. To assess the completeness of our survey, we adopt an empirical approach that involves inserting synthetic sources into our processed data cube. We evaluate their recovery rate by requiring that they be detected via the same process as the real detections in our survey. As described in Section 2.2, and given the number of artificial source injections required to accurately sample the multidimensional parameter space of real H I sources, employing this approach requires a highly automated source-finder that can detect galaxies with high reliability.

All injection/recovery methods are based on prior knowledge of the underlying distributions of source parameters, such as mass, inclination, size, and velocity width. To create a catalog of synthetic sources that accurately samples the multidimensional parameter space spanned by H I sources, we follow an approach based on the one detailed in A. R. Gogate et al. (2020). The synthetic H I sources are simulated using the Groningen Image Processing System (GIPSY; J. M. van der Hulst et al. 1992) package's galmod task, which uses a tilted ring model (D. H. Rogstad et al. 1974) to create the 3D velocity field of a given H I source. For each galaxy, galmod requires the radial H I surface density distribution (in cm -2 ) and the rotational velocity distribution (in km s -1 ) as a function of the radius (in arcseconds) of the galaxy, as well as the velocity dispersion (in km s -1 ), the inclination angle, and the position angle (in degrees). We use the radial surface density distribution profiles as described in J. L. Sérsic (1968) and T. P. K. Martinsson et al. (2016) and the rotational velocity profiles described by M. Persic et al. (1996) and S. Courteau (1997), as applicable to different mass ranges, to generate synthetic H I sources across the wide mass range of our study. We make use of the known local H I size-mass scaling relation (J. Wang et al. 2016) when generating the parameters required by galmod to guarantee that the sources follow this relation. The equation for the amplitude of the rotational velocity profiles (from M. Persic et al. 1996) likewise yields results that are consistent with the Tully-Fisher relation (R. B. Tully & J. R. Fisher 1977). Our final catalog comprises over 150,000 synthetic sources, with randomly chosen inclinations and position angles, covering H I masses in the range 7 log (M H I /M e ) < 10.75, such that the vast majority of sources have masses in the lower half of the range, as required to ensure reliable recovery statistics for the entire mass range (see below). To ensure minimal alterations to the noise properties of the data cube, which can influence the source-finding process, the number of sources in any single injection/recovery trial is limited to a maximum of 500 for sources near the low-mass limit and a maximum of 100 for sources near the high-mass limit.

When SoFiA is run on a given cube that includes both real and synthetic sources, it delivers a list of all spectral line detections and-for each detection-a 3D "cubelet" that includes generous spatial and spectral buffers around the pixels in which emission is seen. 29 We calculate the total H I mass of each detection as

where D L is the cosmological luminosity distance to the source, and S is the integrated H I flux density (M. Meyer et al. 2017).

The integrated H I flux density is calculated using zeroth moment maps that are created by first smoothing individual source cubelets to a circular beam of 20″ × 20″. Each smoothed cubelet is clipped at a 3σ threshold (σ is estimated by measuring the rms noise over an emission-free region) to create a mask that encompasses diffuse low-column-density emission.

We remove isolated regions corresponding to noise peaks from the mask by discarding all regions whose areas are smaller than that of the smoothed beam. We then apply the resulting mask to the original-resolution cubelet to generate the zeroth moment map and calculate the integrated H I flux density. False-positive detections can become a major source of error in recovery rate calculations. Given our perfect knowledge of where the synthetic sources are, we can easily identify any false positives that SoFiA "recovers"; however, it is not possible to similarly distinguish false from real sources in our full observed sample. In order to eliminate the complications that would be caused by false-positive detections in our catalog of real sources, we need to work only with a high-purity subset of SoFiA detections in both the observed and synthetic data. Recovery tests on the synthetic sources have shown that using SoFiA's reliability threshold of 99% creates a high-purity subsample (purity > 99%). 30 All of the real sources in the selected subsample are detected in every injection/recovery trial run, which confirms that the injected population in each trial does not alter source-finding accuracy. The outputs of our 500 trials are then combined to create a reference recovery catalog in which the fate of each injected source and its properties are recorded. We split the total mass range (7 log (M H I /M e ) 11) into eight bins in order to investigate the recovery details for each bin. A traditional 30 Because no analogs of multiwavelength counterparts exist for synthetic sources, in contrast to real sources, we cannot incorporate the use of multiwavelength catalogs into the source-finding process for real sources. Without using multiwavelength data, it becomes much more challenging to determine how many source candidates identified by SoFiA (e.g., above a lower reliability threshold of 95%) are false detections. Therefore, we need to establish a reliability threshold using injection recovery tests alone that delivers high-purity samples (e.g., the high-purity sample from the 99% reliability threshold), which do not require cross-validation against multiwavelength data. 31 The sources are spectrally smoothed in galmod with a Gaussian kernel whose FWHM is twice the channel separation. 32 Accurately capturing the effects of large-scale structure (LSS) at the point of source injection would (ideally) entail simultaneously solving for both the underlying LSS and the HIMF using our observed mass and spatial distributions. Unfortunately, our sample size is insufficient to support such a joint analysis. Instead, we have assessed source recovery statistics in the presence of LSS analogous to that observed in our real source catalog. These tests were focused on sources with

, which represent 75% of our observed sources and sample volumes where a reasonable estimate of LSS can be derived from H I data. Our tests confirm that the presence of LSS in our data along the line of sight does not significantly affect source recovery statistics as a function of H I mass, although other surveys might be affected differently. 33 The spatial position of a given synthetic source is determined in (r, θ) coordinates with respect to the pointing center. Radius r is sampled according to p(r)dr = rdr, and θ is sampled from a uniform distribution between 0 and 2π. Each (r i , θ i , z) tuple is cross-checked with existing source positions to avoid overlap. 34 99% of the injected sources are recovered within a sphere with a 3.5 spaxel radius centered on the injected position.

approach here would calculate the ratio of recovered-toinjected sources in each bin and correct the number of real detections by the recovery factor to determine the underlying distribution of galaxies. However, our simulations show that even with a recovery fraction of close to 100%, a considerable number of galaxies injected in a given mass bin are recovered in a different bin, with a shift into a higher-mass bin much more likely than the alternative. Given that there are significantly more low-mass objects than higher-mass objects in the real Universe, the recovery of sources in higher-mass bins can introduce a bias in the inferred HIMF that needs to be corrected.

In order to keep track of the migration of recovered sources across mass bins, we create a correction matrix-a 2D array of recovery factors (instead of a correction vector, i.e., a 1D array of recovery factors)-in which the fates of sources injected in each mass bin are recorded as a function of their recovered mass. This "recovery matrix" can be thought of as a function over an M inj(ected) versus M rec(overed) domain, where the function at a given (M inj , M rec ) point represents the fraction of sources with an injected mass in the M inj bin that have been recovered with a mass in the M rec bin. In an ideal case, the recovery fraction matrix would be the identity matrix, signifying that each source has been recovered in the same mass bin in which it was injected. Our tests show that the recovery fraction matrix for our data cube is a band-diagonal matrix with significant nonzero values on the main diagonal and the adjacent diagonals, such that sources are most likely to be recovered in their original injection bin and much more likely to be recovered in the next-higher mass bin than in the next-lower mass bin (Figure 2). The greater likelihood of recovering sources in the adjacent higher-mass bin results from the combined effect of several factors, notably the effects of the source-finding process and image-plane continuum subtraction. Source finding relies heavily on threshold cuts, such that a source near the detection limit might be selected as a detection if it coincides with a positive noise region, but ignored if it coincides with a negative noise region. 35 Estimation of the underlying continuum during pixel-based continuum subtraction in turn could be biased toward underestimating the continuum level for some sources depending on their parameters (such as inclination), leading to under-subtraction of continuum and overestimated line fluxes for those sources. The combination of these factors, along with the effects of other steps in our data processing, leads to an asymmetric bias toward recovery at higher masses.

In order to make sure that we have an accurate estimate of each matrix element, we have injected enough sources to ensure that the fractional uncertainties on the diagonal and below-diagonal matrix elements are <10%. Figure 3 shows the percentage uncertainties on the diagonal matrix elements as a function of the number of injected sources. To calculate the uncertainties on the matrix elements for the recovery of synthetic H I sources, we employ a method inspired by the bootstrapping technique. For each mass bin in the analysis, we start with a large pool of synthetic sources-specifically, at least 15% more sources than the maximum number used in any calculation. We perform multiple random draws of synthetic source samples from this pool, calculate matrix elements for each draw, and then estimate the uncertainties. We begin with a small sample size, n = 100, and conduct 4096 random draws of this size from the full synthetic catalog, calculating matrix elements for each draw. The uncertainty for this sample size is then derived from the results by calculating the range that includes 68% of the data around the mean. This procedure is repeated with increasing sample sizes (n = 200, 400, L ), until we reach a size that is within 15% of the size of the total pool (e.g., 10,000 for a pool of 11,500 sources). Each draw is made without replacement to avoid biases that might arise from reusing the same sources within a single estimation round. This stepwise approach allows us to observe how uncertainties in For the lower-mass bins, the off-diagonal elements are large relative to the ondiagonal elements. If a source is not recovered in the injected mass bin, it is more likely to be recovered with a higher mass than with a lower mass.

Figure 3. 1σ fractional uncertainty (as a percent) on the diagonal elements of the recovery fraction matrix shown in Figure 2 as a function of the number of injections N. The fractional uncertainty is determined by calculating the recovery matrix elements for 4096 randomly chosen samples of N sources in a given M inj bin from a larger pool of injected sources in that bin. The final 1σ fractional uncertainties on all of the diagonal (shown here) and below-diagonal matrix elements are <10%.

the recovery matrix elements vary with changes in the sample size (Figure 3), ensuring that the calculated uncertainties are representative of the true variability expected in different sampling scenarios.

Numerous studies (e.g., M. A. Zwaan et al. 2005; M. G. Jones et al. 2018) have demonstrated that the HIMF can be effectively characterized by a P. Schechter (1976) function, defined as

, 3

whose three free parameters are the normalization constant f * , the "knee" mass M * , and the low-mass slope α. We adopt a forwardmodeling approach to estimate the Schechter function parameters that are consistent with our observed data. For each set of parameters (f * , M * , α), we calculate the intrinsic number of H I sources for each bin based on the volume of the survey and the bin widths. The intrinsic numbers for all bins (for a given set of Schechter function parameters) are then multiplied by the rows of the recovery matrix and summed across columns to determine the expected numbers of observed sources in all bins. The likelihood of observing the actual numbers of sources is then calculated based on these expected counts for each set of Schechter parameters during the fitting and Markov Chain Monte Carlo (MCMC) sampling process. By directly estimating the Schechter function parameters using the observed set of detections, we avoid the need to create an unbiased (corrected, binned) histogram of detections for fitting with a Schechter function.

Our forward-modeling approach allows us to convert the numbers of galaxies in the different M inj bins predicted for a given set of Schechter function parameters into numbers of recovered galaxies over the M rec bins. We calculate the total number of recovered galaxies in each M rec bin by summing the number of recovered galaxies from all of the M inj bins. Our tests confirm that the inferred low-mass slope of the Schechter function can be biased if we fail to account for the recovery of sources with higher masses than the bins in which they were injected. Given that nontrivial numbers of sources are recovered with higher-than-injected masses, we have chosen to exclude the lowest M rec bin from the forward-modeling process, as we cannot properly estimate the number of sources with

< that have been recovered in that bin. At the high-mass end, our analysis includes the

bin, in which we actually detect no real H I sources, in recognition of the fact that a nonnegligible fraction of sources are recovered with higher masses and in order to determine the knee mass more accurately. Including the highest-mass bin in the forwardmodeling process constrains the number of sources in the second-highest mass bin to be consistent with the lack of detections in the highest mass bin. We note that the predicted number of sources in a given bin is the mean of a Poisson distribution describing the source count for that bin. The Poisson distribution P(μ) for mean values μ > 7 can be approximated well with a Gaussian distribution ( ) N , mm . However, the discrepancy between the two distributions for smaller values of μ introduces a bias in the forward-modeling process by differently weighting the estimated likelihoods for bins with few detections. We choose the Poisson likelihood, as our investigations reveal that choosing the Gaussian or meanstandard-error likelihood instead of the Poisson likelihood significantly underestimates the uncertainties associated with the inferred Schechter function parameters and introduces a bias toward steeper low-mass slopes.foot_6

It is necessary to maintain consistency between the observed HIMF and the hypothetical HIMF used in completeness calculations. For a given mass bin i, the completeness fraction C i represents the fraction of all sources in that bin that can be detected within the survey volume. In a traditional injection/ recovery method, for bin i, the completeness fraction C i is defined as the ratio of the number of recovered sources R i to the number of injected sources I i in that bin, i.e., C i = R i /I i . To account for cross-bin contamination, we can write

where f j→i represents the fraction of sources that are injected in bin j and recovered in bin i. Isolating the term for galaxies both injected and recovered in bin i results in

showing that when some fraction of sources recovered in bin i are from bin j (for j ≠ i), the completeness fraction for bin i depends on the ratio of the number of injected galaxies in bin j to bin i, i.e., I j /I i . In a traditional approach, the number of injections in each bin is proportional to the predicted count from the HIMF to preserve relative ratios and enhance statistical reliability. For example, if the HIMF predicts N sources in a given bin, then the number of injected sources is some multiple of N. Given that any HIMF varies significantly-by >2 orders of magnitude-across the mass range of a sample like ours, for every source in the highest-mass bin, about 300 sources need to be injected into the lowest-mass bin. This requirement poses a practical challenge for any traditional approach; for example, to obtain reliable statistics for our highest-mass bin if it were to contain only 100 sources, over 30,000 sources would need to be injected into the lowestmass bin alone to maintain ratios consistent with the HIMF. Using a recovery matrix, however, eliminates this constraint by incorporating cross-bin contamination directly into the matrix, allowing the number of injections in each bin to be chosen independently of the other bins. We make use of this flexibility when determining the number of injections by requiring a 10% threshold on the uncertainties associated with the diagonal and below-diagonal matrix elements.

In order to calculate the best-fit Schechter function parameters and their associated uncertainties, we first maximize the Poisson likelihood for our observed data to calculate a set of best-fit parameters, which are then used in an MCMC sampling process to estimate the uncertainties on those parameters. As discussed above, we inject enough sources in each M inj bin (independent from the other M inj bins) to estimate the diagonal and below-diagonal elements of the recovery fraction matrix with <10% uncertainty. However, the estimated matrix elements and their uncertainties depend in part on the slope of the mass distribution of the injected sources within a given bin. A discrepancy between the observed HIMF and the assumed slope of the mass distribution of sources in a given bin, if present, will result in a biased estimation of the Schechter function parameters. To mitigate this effect, we determine the best-fit Schechter function parameters using two iterations. In the first iteration, an initial recovery matrix is calculated using the ALFALFA α.100 Schechter function parameters (M. G. Jones et al. 2018). This matrix is used to maximize the Poisson likelihood for our data to find an initial set of best-fit Schechter function parameters. In the second iteration, we update the slope of the mass distribution of sources in each bin to be the slope recovered from the first iteration and recalculate the matrix elements, which are then used to find a revised set of best-fit Schechter function parameters. The Schechter function parameters from the second iteration change by <1% compared to the results from the first iteration, eliminating the need for any further iterations. 37We use the PyMCfoot_8 implementation of the MCMC method to estimate the uncertainties in the best-fit Schechter function parameters and their covariance by sampling an appropriate posterior probability distribution. The prior distributions for the Schechter function parameters are log 10 (f * ) uniform in [-4, -2], log 10 (M * ) uniform in [9, 10.5], and α uniform in [-1.8, -0.5]. In the MCMC sampling process, we create 512 different realizations of the recovery matrix using the best-fit Schechter function parameters, in order to account for the uncertainties in the matrix elements, and perform an MCMC sampling with eight chains and 32,768 steps (after burn-in) for each of these 512 matrices. The posteriors from the 512 realizations are then combined to create a full posterior, for estimating the uncertainties on the best-fit Schechter function parameters in a way that includes the uncertainties in the matrix elements.

Taken together, the elements of the recovery matrix method for HIMF determination-starting with source injection and recovery as laid out in Section 3.1.1, and ending with the estimation of Schechter function parameters as described immediately above-could potentially be subjected to a full end-to-end validation test. In such a test, ideally, a single realization of a mock galaxy population with a known HIMF could be simulated into a source-free version of the LADUMA DR1 data cube, and the Schechter function parameters recovered from that mock population could be compared to those of the input HIMF. Unfortunately, as discussed below in Section 5.4 and Appendix A, the appropriate (Poisson) uncertainties in the Schechter function parameters for a sample of the size being analyzed in this paper are very large, such that the parameters recovered for a single realization are likely to be very different from the input parameters a priori. As a result, (in)consistency between input and output HIMF parameters for a single mock data cube cannot be used as a test of the recovery matrix method per se. End-to-end validation of the method could in principle be achieved with an ensemble of mock data cubes; however, creating an ensemble of cubes that is sufficiently large to test the recovery matrix method (beyond the ability of Poisson uncertainties to compromise the test) would be computationally prohibitive and is beyond the scope of this paper.

In addition to the forward-modeling approach discussed above, we have calculated the HIMF using the modified maximum likelihood (MML) method described in D. Obreschkow et al. (2018). The MML method, which assumes Poisson statistics, recovers the HIMF without any binning while dealing with mass uncertainties and the selection biases present in the data. This method uses as an input a selection function that describes the recovery rate of sources as a function of distance and mass, which is averaged over extra variables such as width (D. Obreschkow et al. 2018). We estimate the selection function for our data cube by creating 2D histograms of the injected and recovered sources on a distance-mass grid. We have chosen 32 distance bins and 16 mass bins spanning the full distance and mass ranges of our data. The selection function is defined using a linear interpolator (in lieu of an analytical expression) on the 2D mass-distance gridded data, where the gridded values are calculated as ratios of the recovered source histogram and the injected source histogram. We use the dftools package developed in R and described in detail in D. Obreschkow et al. (2018) to calculate the corrected number density of H I sources in each bin. We then fit a Schechter function to the corrected number densities and estimate the uncertainties in the parameters of the Schechter function using MCMC sampling.

The MML method offers the flexibility to calculate the Schechter function parameters in several different ways, such as providing the maximum volume V max in which each source can be detected, or alternatively, providing a selection function. 39 While using the maximum volume option results in a reasonable fit to our data, we chose to generate a selection function based on our extensive injection/recovery tests for a fair comparison with the recovery matrix method, as described in Section 3.2.foot_10 However, this approach results in convergence issues with the MML method, resulting in an unreliable fit. To work around this problem, we use the built-in functionality of the dftools package to extract the corrected HIMF values for our mass bins and subsequently fit a Schechter function to these values. This process closely aligns with previous HIMF studies, where the corrected HIMF is calculated from binned data and a Schechter function is then fit to those values, making this discussion relevant to previous works that follow a similar analytical framework.

In Figure 4, we compare the Schechter functions preferred by the recovery matrix method and the MML method for the LADUMA DR1 data, along with the actual detected number of galaxies in each bin within our high-purity sample. The best-fit parameters for both methods are determined using a search on a variable resolution grid with a high density of points around the best-fit region. Figure 5 shows the estimated Schechter function parameters and their uncertainties plotted on the marginalized 1D and 2D posterior probability distributions for the recovery matrix method. The associated uncertainties in the parameters are determined directly from the 3D sampled posterior, as none of the three Schechter function parameters are considered nuisance parameters. 41 Determining the appropriate uncertainties for the MML result is complicated by the fact that the MML method does not provide uncertainties for any bin with zero detections. If we circumvent this limitation by excluding the zero-detection highest-mass bin, we obtain only a weak constraint on the knee mass, with a 1σ upper limit reaching the upper limit ( detections in that bin (only modestly higher than the actual value of 0, and encompassing both 0 and 1 within its 1σ Poisson uncertainty). This adjustment to the highest-mass bin leads to 41 The 4D volume enclosed inside the 1σ limit of a 3D Gaussian distribution is 24.91% of the total volume. Estimating the parameters and their uncertainties from the marginalized posterior distributions does not lead to the same results as when parameters are estimated from the 3D distribution (see Figure 5).

more reliable estimates of the knee mass and its associated uncertainties. Given that the output of the MML method is the corrected number density and associated uncertainty for each bin, we need to use a Gaussian likelihood for fitting and MCMC sampling. Therefore, the MML results are slightly biased by the effects of bins with low numbers of counts (as explained in Appendix A) compared to the results from the recovery matrix method. The results from the MML method, which are consistent with our recovery matrix method results but less so with the results of previous surveys, are shown in Figure 6.

Comparing our results with those of previous surveys is complicated by the fact that different authors have used different values of H 0 (with which f * , M * , and Ω H I scale straightforwardly), different approaches to defining cosmological volume (which are not always stated but will affect f * and Ω H I ), and different choices of likelihood that translate to different "1σ" uncertainties. To provide as consistent comparisons with previous results as we can, in Table 2 we have scaled the values of f * , M * , and Ω H I reported by the HIPASS, AUDS, and MIGHTEE teams for our choice of H 0 (ALFALFA uses the same value of H 0 , so requires no rescaling), and we provide indications in the table notes of how we might expect results to change further on the basis of volume calculations. We also present the best-fit Schechter function parameters for the recovery matrix method and the MML method for LADUMA, along with a "free fit" version of the MML results that illustrates the impact of ignoring the zero-detection highest-mass bin. Figure 7 shows that the estimated Schechter function parameters from the recovery matrix and MML methods are in good agreement with each other. Relative to previous HIMF measurements from the literature, while the individual Schechter function parameters that we recover with the MML method are not consistent with all previous results (e.g., the AUDS measurement for α), those we recover for LADUMA using our (preferred) recovery matrix method agree with those for all previous surveys within the associated 1σ uncertainties when added in quadrature. The situation for the HIMF as a whole is not as clear: due to covariance between the Schechter function parameters, some of the literature HIMFs fall outside the 1σ uncertainty swath shown in Figure 4 across a wide range of masses. A full assessment of the consistency of our derived HIMF with the results of previous surveys, which would require detailed knowledge of the covariances among those surveys' respective Schechter function parameters and perhaps recalculation of their uncertainties using a Poisson likelihood (which then might or might not overlap with the LADUMA uncertainty range), is beyond the scope of this paper.

A complete inventory of H I in the local Universe would include neutral hydrogen that lies outside of galaxies, both within and beyond the cosmic web. In this paper, we can estimate the contribution of galaxies to Ω H I based on the comoving H I mass density (ρ H I ) that is calculated by integrating the Schechter function. The density ρ H I is computed analytically (M. Meyer et al. 2017) as ( ) ( ) M 2 , 6 H I ra f =G + * * where is the Euler gamma function and f * , M * , and α are the Schechter function parameters. The contribution of galaxies to Ω H I is then calculated as ( ) G H 8 3 , 7 H I 0 2 H I p r W= where G is the gravitational constant, and H 0 is the Hubble constant. We measure a galactic Ω H I = 3.09 10 0.49 0.58 4 -+ -with the recovery matrix method and a galactic Ω H I = 3.48 0.93 0.26 -+ 10 4 -with the MML method. Uncertainties in these measurements encompass 68% of Ω H I values (around the best-fit value) calculated from 4096 sets of Schechter function parameters randomly drawn from the 3D sampled posterior. Our results are in agreement with the results for all of the previous surveys except for MIGHTEE within 1σ uncertainties added in quadrature.

The HIMF we have measured from the LADUMA data is subject to some of the same caveats as mass functions measured from previous surveys, in ways that are instructive about the limitations that apply to any HIMF measurement. One significant caveat is that the detectability of a source depends on its H I velocity width as well as its H I mass; our recovery matrix approach deals with this complication by marginalizing over the expected distribution of disk inclinations (see, e.g., D. Obreschkow et al. 2018). Perhaps the most important caveat is that for a fluxlimited H I sample, the sources in the lowest-mass bin are typically confined to a much smaller (nearby) volume than the total survey volume (M. A. Zwaan et al. 1997; see especially their Figure 4). While previous authors have emphasized the importance of the bias of a high-mass galaxy population relative to the underlying dark matter distribution as a contributor to cosmic variance (e.g., B. P. Moster et al. 2011), the limited volume in which low-mass galaxies can be detected for a flux-limited sample can be just as important. Among the three Schechter function parameters, M * will be least sensitive to "available" volume and most sensitive to population bias, as it is predominantly determined by the high-mass bins. In contrast, and as pointed out by M. G. Jones et al. (2018), α is expected to be most sensitive to available volume (albeit least sensitive to population bias), as it is predominantly determined by the low-mass bins. f * is expected to fall somewhere in between in terms of susceptibility to cosmic variance, as all bins contribute to the overall normalization within the limits of their associated uncertainties. For our data, the accessible volume associated with the lowest-mass bin, in which we recover any sources in our highpurity sample (specifically, two sources), is only about 2% of the total survey volume. We find that including this bin in our calculations for the best-fit Schechter function parameters would result in a significantly steeper value for α. Our decision to exclude this bin from our calculations reduces the impact of the extremely small volume of the lowest mass bin. Appendix B uses the ALFALFA α. 100 galaxy catalog to demonstrate that cosmic variance increases the uncertainty in the number of high-mass

) galaxies in the LADUMA volume compared to Poisson uncertainties alone, but the accessible volume subtleties discussed above preclude any extension of this analysis to individual Schechter function parameters.

We note that LADUMA does derive some benefit for HIMF derivation from its greater depth compared to previous surveys: as Figure 8 illustrates, the volume that is available to LADUMA for constraining the HIMF at

is only ~1 order of magnitude smaller than those that were available to HIPASS and ALFALFA, in contrast to the ~4 orders of magnitude advantage those surveys have over LADUMA in terms of available volume for the highest

) masses. From this point of view, the fact that LADUMA recovers a faint-end slope α flatter than (but consistent within the errors with) ALFALFA-and indeed matching the ALFALFA measurements of α in certain of that survey's subvolumes (M. G. Jones et al. 2018)-is not an unexpected result. At higher masses, LADUMA benefits in a different way from its greater depth, namely superior resilience against the systematic loss of high-inclination, high-mass sources at large distances within its observed sample, as demonstrated in Figure 9.

The behavior of the observed LADUMA HIMF (see Figure 6), particularly the rolloff just below M * followed by an upturn at lower masses, is reminiscent of features seen in optical/near-IR luminosity functions (LFs). These LFs can be effectively decomposed by galaxy type or color into multiple LFs, each represented by its own Schechter function (A. Sandage et al. 1985;J. Loveday et al. 2012;S. P. Driver et al. 2022). Previous studies (e.g., K. Said et al. 2019;M. G. Jones et al. 2020) have shown a dependency of the Schechter function parameters on environment, suggesting that a more complex Notes. z max shows the upper limit of the redshift range for each survey. For LADUMA, the disfavored MML "free fit" shows the best-fit values obtained by excluding the highest-mass bin, which has no detections, from the fitting process. We calculate the comoving volumes using the area and redshift range reported for each survey in the corresponding publication, scaled as needed to our cosmology. All other parameters are rescaled as needed for H 0 = 70 km s -1 Mpc -1 . References in the table are as follows: J18 = M. G. Jones et al. (2018), Z05 = M. A. Zwaan et al. (2005), X21 = H. Xi et al. (2021), and P23 = A. A. Ponomareva et al. (2023). a We list here the comoving volume probed by the ALFALFA α.100 sample for our assumed cosmology, which is slightly smaller than the 6.5 × 10 6 Mpc 3 actually used in the HIMF calculations of M. G.

Jones et al. (2018; see also M. Jones 2024, private communication). Use of our smaller volume would modestly increase f * and Ω H I for ALFALFA relative to the values reported in this table. b The reported ALFALFA Ω H I is corrected for H I self-absorption; the value without this correction is 3.50 0.51 0.51 -+

. c We list here the comoving volume probed by the MIGHTEE analysis of the COSMOS and XMM-LSS fields for our assumed cosmology, which is (a) larger than the 7 × 10 3 Mpc 3 "cosmological volume" reported in A. A. Ponomareva et al. (2023) that includes an unnecessary extra factor of h -3 , but (b) only slightly different from the volume actually used in that paper's V eff analysis (A. Ponomareva 2024, private communication). Use of our slightly different survey volume in the context of a V eff analysis would modestly change f * and Ω H I for MIGHTEE relative to the values reported in this table. LADUMA's f * and M * agree with the results for all of the previous surveys within 1σ uncertainties. LADUMA's low-mass slope is shallower than those measured in all previous surveys, although it is in agreement within the 1σ uncertainty with ALFALFA and MIGHTEE results. The values for these parameters are reported in Table 2, along with the literature references for previous surveys. The best-fit parameters from previous surveys have been rescaled to H 0 = 70 km s -1 Mpc -1 . functional form (e.g., different Schechter functions for highdensity versus low-density environments) might be required to describe the overall HIMF accurately. Our forward-modeling approach is designed to accommodate various forms of mass functions beyond the traditional Schechter form, providing a flexible framework for analyzing the HIMF in greater detail for larger data sets. In addition to more complex functional forms, our approach can easily incorporate nonparametric models by directly sampling the intrinsic numbers of galaxies in different bins based on different priors, predicting corrected number densities for all bins. However, our current data set does not permit the fitting of models with more degrees of freedom, as the relatively small number of sources available cannot adequately constrain a more complex model.

This Section and the subsequent Section highlight two challenges in HIMF determination, illustrating the merits and limitations of the recovery matrix approach compared to traditional methods. The completeness of a survey for sources of a certain mass can be assessed through two primary methods. The effective volume (V eff ) method estimates the survey's completeness for sources of mass M by evaluating the effective volume within which these sources can be observed.foot_11 Injection/recovery methods offer an alternative approach by introducing synthetic sources into the data cube and tracking their detection or nondetection using the same processes applied to real sources. These approaches can be encapsulated in expressions for the measured number density of sources of mass M, f(M),

where n o (M) represents the number of sources detected in a bin of width d M log centered at M, V eff (M) is the effective volume for sources with H I mass M, V s is the total volume of the survey, and C s (M) is the completeness of the survey for sources with H I mass M (estimated via injection/recovery tests). By comparing Equations (8) and (9), it becomes evident that the effective volume is equivalent to the completeness multiplied by the total volume of the survey. The completeness C s can be further as

, where C p (M) is the completeness based on the distribution of source properties, such as inclination, and C V (M) is the fraction of the total volume of the survey available for detection of sources of mass M. As described in Section 5.1 and demonstrated in Figure 8, the accessible volume for sources of mass M can vary drastically, e.g., by orders of magnitude, over the mass range of a survey. The completeness for low-mass sources is largely determined by the fraction of the total volume accessible to those sources and less so by the distribution of their intrinsic properties. The combination of these factors reveals that fluctuations (e.g., in the number of observed sources, which is sensitive to the distribution of the sources) in the small volume available to sources at the low-mass end of the mass range are amplified by the inverse of Figure 8. Volumes accessible for the detection of sources with given M H I by LADUMA, ALFALFA, and HIPASS, calculated based on their respective detection limits. The maximum volume inside which all sources of H I mass M H I are detected cannot be precisely defined, as it depends on parameters such as velocity width. However, we can roughly define an "accessible" volume beyond which a given survey in practice detects vanishingly sources of H I mass M H I ; these are the volumes plotted above. The number density distribution for sources of H I mass M H I is sensitive to the distribution of sources (large-scale structure) inside the corresponding "accessible" volume. While all surveys have a much smaller (by several orders of magnitudes) available volume at lower masses compared to higher masses, that discrepancy is less significant for LADUMA compared to HIPASS and ALFALFA. The fact that the flat part of the LADUMA curve (corresponding to its total survey volume) extends farther left than do the flat parts of the HIPASS and ALFALFA curves is due to LADUMA's greater depth (and in spite of LADUMA's slightly larger redshift range) relative to the other two surveys.

Figure 9. Recovery fraction as a function of source inclination for the higherredshift half of LADUMA's low-z SPW. While the recovery fraction for lowmass sources decreases significantly at high inclinations, the recovery fraction for high-mass sources remains consistently high at all inclinations.

the accessible volume fraction in the completeness correction process. Therefore, it is essential to limit the amplification of these fluctuations, which can be done by imposing a minimum accepted volume fraction (10% for this work) for the low versus the high end of the mass range.

Beyond the amplification by the inverse of the volume fraction, as laid out in detail in Appendix C, the accuracy of the completeness correction in the V eff method is sensitive to the distribution of properties in the observed population. To demonstrate this point, we consider two different typesfoot_12 of sources of mass M: face-on and edge-on. For unresolved detections, it is clear that edge-on sources can only be detected in a smaller volume than otherwise identical face-on sources, since their larger line-of-sight velocity widths translate to lower peak flux densities. Let us represent the ratio of the volume available to face-on sources to the volume available to edge-on source as R f/e . Assuming (unrealistically) that there are intrinsically the same numbers of edge-on and face-on sources in a survey volume, we would expect that for each edge-on source we detect, R f/e face-on sources should be detected. As described in Appendix C, the effective volume calculated for sources of mass M is sensitive to the observed ratio of the number of face-on to the number of edge-on sources. However, for an injection/recovery approach, this sensitivity would not affect the inferred completeness of the survey as long as the underlying distribution of sources is well understood,foot_13 since the completeness is determined through source injection.

Another challenge for HIMF determination involves accurately addressing uncertainties in the fitting process, particularly in bins with few or no detections. The data consist of observed source counts across various mass bins, each representing a count randomly drawn from a Poisson distribution. A traditional approach assumes that the counts are Gaussian-distributed with a Poisson uncertainty-computed as the square root of the number of sources in each bin-and uses a completeness correction factor to calculate the corrected number density of H I sources per bin. These values are then used to estimate the best-fit Schechter function parameters and their uncertainties. As outlined in Appendix A, using a Gaussian approximation for the full Poisson likelihood significantly biases the derived best-fit parameters of the Schechter function and tends to underestimate the 1σ uncertainties, especially in bins with few detections. This underestimate affects both bins where completeness is high but galaxies are few in number, i.e., at high masses, and bins where galaxies are plentiful but completeness is low, i.e., at low masses. 45 The root of this discrepancy lies in the asymmetric and positively skewed nature of the Poisson distribution. For a bin with a small number of detections, assuming Gaussiandistributed counts skews the estimated mean toward the observed count, leading to overfitting of the data. For LADUMA, this bias would manifest as an unrealistically steep low-mass slope, since (as is typical for any survey) the lowest-mass bins contain relatively few detections. In contrast to a traditional approach, the recovery matrix method employs forward modeling to predict the mean of the underlying Poisson distribution for each mass bin for any set of Schechter function parameters and implements a Poisson likelihood for model fitting and MCMC sampling. One significant advantage of forward modeling is its ability to integrate bins with zero detections into the fitting process, enhancing the robustness of the fit. As discussed in Section 4.1, for the MML method, we observe significantly larger uncertainties in the knee mass parameter when we exclude the zero-detection highest-mass bin. However, forward modeling naturally accommodates such cases, as it models uncertainties based on the predicted values and is compatible with the assumption of a Poisson distribution for the observed counts.

In this paper, we present a comprehensive analysis of the HIMF and the contribution of galaxies to Ω H I from a portion of the LADUMA survey covering 0 z 0.088. Using a high-purity sample of H I detections from LADUMA's DR1, we develop a new "recovery matrix" method and benchmark it against a traditional maximum likelihood approach for measuring the HIMF. Our estimates of the

Schechter function parameters are 3.56 10 1.51 1.79 3 f =-+ -Mpc -3 dex -1 , 1.18 0.14 0.14 a =--+ , and  ( ) M M log 10.01 0.17 0.23 = -+

. These values are in agreement with previous measurements of the HIMF over similar redshift ranges, enhancing confidence in the robustness of our results and the resilience of LADUMA's single-pointing geometry against the effects of cosmic variance.

Our methodology uses extensive synthetic source injections to correct for survey incompleteness and, for the first time, includes the effects of the continuum subtraction process. We account for varying sensitivity and completeness across the survey volume and mass range. In particular, our forward-modeling approach proves beneficial in handling bins with few or no detections (using an appropriate Poisson likelihood during the fitting and MCMC sampling process), thereby minimizing systematic effects on the derived Schechter function parameters. By using a Poisson likelihood instead of a Gaussian likelihood with Poisson uncertainties, we avoid overfitting the low-mass and high-mass bins of the Schechter function, which usually have few detections. Our analysis highlights the importance of cultivating a high-purity sample for reliably estimating survey completeness and avoiding biases in the recovery of mass function parameters. The recovery matrix allows for independent determination of the required number of injections per mass bin, ensuring reliable statistics across all bins-overcoming a challenge faced by traditional injection/recovery methods.

Looking forward, the LADUMA survey's future data releases will allow us to refine these measurements and potentially reveal new aspects of H I evolution across a broader redshift range (e.g., M. Hoosain et al. 2025, in preparation). The continued development and application of new methods like the recovery matrix will be crucial in leveraging the full potential of these data to enhance our understanding of galaxies' H I distributions and their role in galaxy evolution. Note. μ i and y i represent the predicted and observed counts, respectively. 46 When errors are calculated based on observed data, the error term in the overall likelihood is a constant-( )

using the terminology of Table 3 -and can be factored out of the likelihood sum.

To provide an indication of how cosmic variance affects our results, we use the ALFALFA α. 100 H I source catalog (M. P. Haynes et al. 2018) to estimate how many high-mass H I galaxies would be detected within volumes comparable to that of LADUMA's low-z SPW cube. We begin by choosing a subsample of the ALFALFA catalog with high completeness, by excluding sources with log(M H I /M e ) < 9.75 and z > 0.05. These cuts minimize the impact of survey incompleteness and result in a sample of over 8500 galaxies. To ensure a fair comparison with the low-z SPW for LADUMA, we define 388 independent "pencil-beam" subvolumes within the ALFALFA volume; each of these is matched to the volume of the LADUMA low-z SPW by selecting a larger solid angle on the sky to compensate for the smaller redshift coverage of ALFALFA. Within each defined subvolume, we count the number of galaxies. We then calculate the standard deviation in galaxy counts across these subvolumes, which is found to be ∼8.5 (relative to a mean of 16). This observed standard deviation exceeds the expected Poisson uncertainty of 16 4 = , providing an indication of the impact of cosmic variance on our results. In our LADUMA sample, the observed number of high-mass sources is 15, closely matching the average number derived from the ALFALFA test. While this exercise suggests that Poisson uncertainties represent only ∼47% of the total uncertainty when cosmic variance is included, the facts that uncertainties are not symmetrically distributed for the Schechter function parameters and that the ALFALFA sample does not provide a sufficiently robust estimate of cosmic variance for low-mass galaxies (due to the smaller available volumes) mean that extending this analysis to individual Schechter function parameters is nontrivial.

Considering the intrinsic number density, f t , of sources with H I mass M, i.e., the value of the HIMF at M H I = M, we can write:

where d M log is the width of a bin centered at M, and n t is the intrinsic number of sources (not necessarily an integer) within a survey's volume, denoted by V s . 47 We can express the observed number density of sources with H I mass M in a survey using a V max (VM) framework as

and using an injection/recovery (IR) framework as

where n o represents the integer number of sources detected in an interval of width d M log centered at M, which can be described as a randomly sampled (integer) value from a Poisson distribution with mean μ = C s n t , i.e., n o ~P(μ = C s n t ); V eff is the effective volume of sources with H I mass M (estimated using the observed sample); and C s is the overall completeness of the survey for all types of sources with H I mass M (estimated via injection/recovery tests). Comparing Equations (C2) and (C3) shows that the V max approach estimates the completeness of a survey by calculating the effective volume (V eff ) using the observed source population. It is clear that any error introduced by the randomness of the observed population affects the numerators in Equations (C2) and (C3) identically. For the denominators, this analogy breaks down. The denominator in Equation (C3) is not affected by errors due to the randomness of the observed population, as V s is fixed and the completeness, C s , is determined via injection/ recovery tests. However, the denominator in Equation (C2)when determined using the observed population-will be affected by errors due to the randomness in the observed sample. It is helpful to track how random Poisson uncertainties affect the observed number density for these two methods. For this exercise, we assume that the uncertainties introduced by the injection/recovery process are negligible (as they can be reduced by increasing the number of injected sources) compared to the random uncertainties.

First, we need to consider the different "types" of galaxies with M H I = M. We assume there are N types of galaxies with normalized abundances of x i , where

x

; one example would be galaxies with the same mass but different inclination angles. Assuming a completeness fraction of C ifoot_16 for each source type, the expected number of sources 〈n o 〉 in the survey volume would be

are, respectively, the intrinsic number density of sources of type i, the effective volume available for the detection of sources of type i, and the completeness of the survey for sources of H I mass M, and 〈n i 〉 is the expected number of sources of type i in the survey volume. Equations (C5) and (C6) are two representations of the expected number of sources: the sum of the expected numbers of sources of type i (C5) and the intrinsic number of sources multiplied by the completeness of the survey (C6). The number of observed sources n o can be expanded as

( ) ( ) n n n n n n n n , C7 o o o i N i i N i i N i i i N i 1 1 1 1 dd d = á ñ+ =S á ñ+S =S á ñ+ =S == ==

where δn o and δn i (which may be positive or negative) are, respectively, the random uncertainty in the total number of observed sources and the random uncertainty in the number of observed sources of type i, and n i is the number of observed sources of type i. 49In an injection/recovery framework, by combining Equations (C3), (C6), and (C7), we have: , the observed number density is equal to the intrinsic number density. Similarly, in a V max framework, we can combine Equations (C2), (C5), and (C7) to get:

which shows that the uncertainty in the measured number density is proportional to the sum of the random uncertainties in the numbers of observed sources with M H I = M of type i divided by the effective volumes for sources of type i. Equation (C13) shows that when

, the measured number density would be equal to the intrinsic number density. The trivial case is when δn i = 0 for all source types, i.e., for all source types, the number of observed sources of type i is equal to expected number of sources of type i, n i = 〈n i 〉. The more complicated scenario will be for all of the nonzero {δn i } in the sum to cancel each other. An interesting case is to compare Equations (C9) and (C13) in the limiting case of a volumelimited sample for all sources of type i. For a volume-limited sample, the effective volume for each source type is equal to the total volume of the survey, i.e., V i = V s , and the completeness for sources of type i is equal to 1, as all such sources can be detected throughout the total volume, i.e., C i = 1. Substitution of these values into Equations (C9) and (C13) results in

making it clear that when all types of sources of H I mass M are detectable throughout the entire volume of the survey, the two approaches are equivalent. We note that in these derivations we have assumed that the effective volume and the completeness of the survey are perfectly determined, underscoring the fact that the difference between these two approaches arises from the methods themselves and are not due to imperfect calculations. In practice, the effective volume and the completeness fraction calculations are affected by errors that contribute to the overall uncertainties in the measured number density distributions. Considering a simplified case of two types of galaxies of equal abundance and identical H I mass M, following an approach similar to Appendix C of D. Obreschkow et al. (2018), we demonstrate how random uncertainties affect the observed number density of sources with M H I = M. Assuming that f t (M) = 2, x 1 = x 2 = 0.5, C 1 = 1/8, C 2 = 1, and V s = 8, we can calculate the expected number of galaxies as

. Distributions of observed number density of sources with M H I = M for the V max and injection/recovery methods. It is clear that the V max method leads to a wider distribution of observed number densities around the true number density of sources with M H I = M. This effect is due to the fact that the V max method uses the observed data set to estimate the completeness of the survey and therefore is affected by uncertainties in the observed population to a larger degree than the injection/recovery method. It is important to note that these differences are inherent in the methods themselves, and this figure assumes perfect knowledge of the completeness of the survey and the effective volume for the (two) different types of galaxies of H I mass M.

( ) 1 8 9. C20 =+= Therefore, we expect to observe one galaxy of type 1 and eight galaxies of type 2, totaling to nine galaxies of M H I = M in the survey volume. We can randomly sample two Poisson distributions with means of μ 1 = 1 and μ 2 = 8 to simulate observations and calculate the observed number density of sources with M H I = M for the two approaches using Equations (C10) and (C8). A histogram of the outcome of 2 13 simulated observations for cases of f t (M) = 2 and f t (M) = 40 is shown in Figure 10, demonstrating that the V max approach leads to a wider distribution of observed number densities around the true number density of sources with M H I = M, due to the fact that the V max approach uses the observed data set to estimate the completeness of the survey and therefore is affected by the uncertainties in the observed population to a larger degree compared to the injection/ recovery approach.

We close by noting that the preceding analysis treats the parameter (M) whose distribution (the HIMF) we are trying to derive in a fundamentally different way from other parameters (e.g., inclination) that determine which "types" different galaxies have, in the sense of affecting the maximum volumes in which they can be detected. This distinction is not absolute: if we were interested (for example) in deriving the H I velocity width function, we could do so treating M as a parameter that determines "type," as it affects the effective available volume for the detection of sources of a given velocity width, and the limitations of the V max approach relative to the injection/ recovery approach would remain the same as we have characterized them above.