Modeling Post-Reionization HI Distributions in Fuzzy Dark Matter Cosmologies Using Conditional Normalizing Flows

Upcoming 21 cm intensity mapping experiments like the Square Kilometer Array (SKA) hold significant potential to constrain the properties of dark matter. In this work, we model neutral hydrogen (HI) distributions using high-resolution hydrodynamical $N$-body simulations of both cold dark matter (CDM) and fuzzy dark matter (FDM) cosmologies in the post-reionization redshift range of $z=3.42-4.94$. We show that the HI abundance decreases in FDM-like cosmologies. Extreme FDM models with $m\sim 10^{-22}$ eV are at odds with a range of measurements. Due to the increased halo bias, the HI bias increases, paralleled by the damped Lyman-$\alpha$ (DLA) bias which we infer from the cross-section of DLAs. The distribution of the latter in extreme FDM models has a high median at the low-mass end, which can be traced to the high column density of cosmic filaments. FDM models exhibit a very similar abundance of DLAs compared to CDM while sub-DLAs are already less abundant. We study the prospects of detecting the brightest HI peaks with SKA1-Low at $z=4.94$, indicating moderate signal-to-noise ratios (SNR) at angular resolution $\theta_A = 2^{\prime}$ with a rapidly declining SNR for lower values of $\theta_{A}$. After training the conditional normalizing flow network HIGlow on 2D HI maps, we interpolate its latent space of axion masses to predict the peak flux for a new, synthetic FDM cosmology, finding good agreement with expectations. This work thus underscores the potential of normalizing flows in capturing complex, non-linear structures within HI maps, offering a versatile tool for conditional sample generation and prediction tasks.


INTRODUCTION
21 cm cosmology studies the redshifted 21 cm-wavelength photons that are emitted in a hyperfine transition of atomic hydrogen atoms (HI).The transition occurs when an electron in the excited spin triplet state flips its spin relative to the proton and falls into the singlet state.While this transition is quantum mechanically 'forbidden' to first order, a very large number of neutral hydrogen atoms in massive clouds render the signal observable.
The hyperfine transition can be used to trace the matter field through a technique called 21 cm intensity mapping (IM) (Chang et al. 2008;Wyithe & Loeb 2008;Santos et al. 2015;Villaescusa-Navarro et al. 2015).It is possible since the spin transition in HI is optically thin (Furlanetto et al. 2006), meaning that 21 cm photons are unlikely to be absorbed/emitted more than once as they travel to our telescopes.Therefore, full 3D mappings of the 21 cm line become possible, with the redshift of the signal providing information about the line of sight distance.Unlike galaxy redshift surveys, there is no need to resolve individual galaxies, which is often required only to determine their redshift.By mapping the unresolved emission of all HI at each frequency of observation and using the ⋆ E-mail: td448@cam.ac.uk observed frequency-redshift relation ν = ν 21 /(1 + z), we can directly map to the corresponding redshift.This is further facilitated by the ease at which spectral resolution can be obtained in radio astronomy.
Ongoing and upcoming 21 cm IM experiments such as HIRAX (Newburgh et al. 2016), PUMA (Bandura et al. 2019), CHIME (Collaboration et al. 2022) and SKA (Combes 2021) are set to yield a wealth of insights into cosmology and astrophysics.With minimal angular resolution requirements and expansive fields of view, these instruments can efficiently map vast cosmological volumes.The 21 cm signal in the post-reionization era (and also during reionization) holds great promise for studying alternative dark matter (DM) scenarios.The 21 cm signal is not only sensitive to DM decays, annihilation processes (Liu & Slatyer 2018;List et al. 2020) or interactions between DM and standard model particles (Barkana 2018;Muñoz & Loeb 2018).It can also probe DM models that result in a suppression of the small-scale matter power spectrum (Sitwell et al. 2014;Schneider 2018;Lidz & Hui 2018;Nebrin et al. 2019;Boyarsky et al. 2019;Jones et al. 2021;Bauer et al. 2021), which are the focus of this paper.
DM models with a small-scale suppression are of interest since the standard model of cosmology and DM -the ΛCDM modelfaces some challenges including the lack of direct detection of DM particles and small-scale observations that are at odds with ΛCDM predictions (Weinberg et al. 2015;Bullock & Boylan-Kolchin 2017).While attempts at unified baryonic solutions to the ΛCDM smallscale problems exist (as reviewed in Del Popolo & Le Delliou 2017), here we investigate fuzzy dark matter (Hu et al. 2000, FDM), which is made up of ultralight bosonic particles of mass ∼ 10 −22 eV called axions, corresponding to de Broglie wavelengths of λ dB ∼ 1 kpc.The motivations behind considering multiple species of light axions is vast, encompassing well-established predictions of string/Mtheory (Arvanitaki et al. 2010;Demirtas et al. 2018) as well as various field theory extensions of the standard model (Peccei & Quinn 1977;Kim & Marsh 2016).One of the distinctive attributes of FDM is its largely redshift-insensitive comoving de Broglie wavelength λ db,c (Khlopov et al. 1985) which scales as λ db,c ∼ (1 + z) 1/4 m −1/2 .This unique property simultaneously addresses challenges associated with small-scale structure and constrains the central density of collapsed halos (Schive et al. 2014), offering a natural resolution to some of the limitations faced by ΛCDM.Notably, FDM introduces solitonic cores within halo centers, which, in models incorporating axion self-interactions, can undergo a phase transition from dilute to denser states (Mocz et al. 2023).
In our investigation, we run N-body hydrodynamical simulations of CDM and FDM cosmologies (as detailed in Section 2.1), exploring a range of axion masses spanning from m = 10 −22 eV to 2×10 −21 eV.Note that in this mass range, FDM is effectively ruled out as comprising 100% of the DM content (as discussed in Dome et al. 2023a).However, in this study, the FDM-like modeling approach will allow us to elucidate trends associated with the axion mass m across a variety of observables.These trends are expected to stay robust even if FDM constitutes only a subcomponent of the DM.
In order to extract the maximum information from IM surveys, it is critical to reliably model the spatial distribution of HI.In our modeling approach, we largely follow the cell-based HI post-processing techniques from Villaescusa-Navarro et al. (2014, 2018); Carucci et al. (2015); Carucci (2018); Diemer et al. (2018Diemer et al. ( , 2019)).We also harness the potential of machine learning (ML) and use normalizing flows, a generative ML model introduced by Agnelli et al. (2010), to capture intricate structures resulting from non-linear physics, facilitating the conditional generation of HI maps with varying axion masses.Specifically, we use a modified version (see Sec. A) of the conditional normalizing flow framework HIGlow (Friedman & Hassan 2022) and showcase the efficiency of the model in sample generation.We assess the generated HI maps against external simulation data, utilizing metrics like the HI mass probability density function and HI power spectrum, validating the prowess of the model across a broad spectrum of mass scales and spatial dimensions.As proof of concept, we demonstrate the ability of HIGlow to interpolate the latent space of axion masses by predicting the peak flux in a synthetic cosmology, see Sec. 4. This affirms its efficacy in characterizing HI distributions for forthcoming parameter forecasting endeavors.
The organization of the paper is as follows: In Sec.2.1, we describe our CDM and FDM hydrodynamical simulations.We summarize our modeling of hydrogen phases and its link to 21 cm physics in Secs.2.2 and 2.3, respectively.In Sec.3.1, we evaluate the trained model.The post-reionization HI abundance, HI and 21 cm clustering (including the HI bias) and the HI column density distribution are investigated in Secs. 3.2,3.3 and 3.4,respectively.For the first time, we analyze the cross-section of damped Lyman-α absorbers (DLAs) in FDM-like cosmologies (in Sec.3.5) and the imaging prospects of SKA-Low using mock radio maps, augmented by the trained normalizing flow (in Sec. 4).We conclude in Sec. 5. HIGlow implementation choices are summarized in App. A.

CDM and FDM-like Simulations
We focus on axions generated via vacuum realignment assuming gravitational interactions do not re-thermalize axions.Simulating such FDM models is significantly more challenging than for CDM, as in the densest regions the wavelike matter oscillations can attain high frequencies ω ∝ m −1 λ −2 dB , requiring very fine temporal resolution even for moderate spatial resolution (Mocz et al. 2019;May & Springel 2021).To bypass the challenges of FDM simulations, here we employ FDM-like modeling described in Dome et al. (2023a).In short, we impose a cutoff in the primordial power spectrum but evolve the system using only the CDM dynamics (also see Ni et al. 2019), providing a useful approximation for FDM and other DM scenarios incorporating a small-scale cutoff such as warm dark matter (WDM, Paduroiu 2022).Following Dome et al. (2023a), we call this proxy classical FDM (cFDM).For scales corresponding to wavenumbers k ∼ 0.16 − 80 h/Mpc, which we explore here, the dynamical manifestation of FDM -the quantum pressure (fluid formulation, Madelung 1927) -has only small impact on the growth of DM fluctuations.Specifically, the absolute fractional difference between growth rates in FDM vs cFDM is less than 5% for particle masses around m ∼ 10 −22 eV and halo mass scales around M ∼ 4×10 9 M ⊙ /h (see e.g.Corasaniti et al. 2017).As opposed to a superfluid, cFDM approximates FDM as a classical collisionless fluid, governed by the Vlasov-Poisson system of equations, but with FDM initial conditions.The exponential-like small-scale suppression in the primordial power spectrum, which we will often refer to as a cutoff, is modeled using the Boltzmann solver AxionCamb (Hložek et al. 2015).
We employ the IllustrisTNG galaxy formation module and run the simulations with the state-of-the-art code Arepo described by Springel (2010) and Weinberger et al. (2020).The hydrodynamical equations are solved on a moving Voronoi mesh using a finite volume technique.Various astrophysical processes such as metal-line cooling, star formation and feedback remain unresolved in the simulations and are approximated by subgrid models (Pillepich et al. 2018).Gas above a density threshold of n H ∼ 0.1 cm −3 spawns star particles stochastically following the empirical Kennicutt-Schmidt relation and assuming a Chabrier (2003) initial mass function (IMF).
In our suite, we use cosmological volumes with a box side length of L box = 40 h −1 Mpc.The box size balances the competing demands of high resolution of HI distributions and large volume (to obtain accurate statistical distributions).We vary the DM resolutions across N = 256 3 , 512 3 , and 1024 3 .For the bulk of this study, our primary focus is on high-resolution runs at 1024 3 while we address resolution effects in our discussions and commentary throughout the text.In such full hydrodynamical runs, the imprints of cFDM are entangled with resolution effects that are due to baryonic physics (e.g.Vogelsberger et al. 2013;Chua et al. 2019).Similar to Villaescusa-Navarro et al. (2018), our findings related to HI distributions are not converged against resolution.Specifically, the HI and 21 cm power spectrum, HI abundance, HI column density distributions and DLA cross-sections exhibit resolution-dependent behavior.However, this does not undermine the validity of our results, as our work is conducted at the effective resolution level of IllustrisTNG100, on which the model parameters have been calibrated to accurately reproduce galaxy properties.Whilst our box size is smaller than that of IllustrisTNG100, our effective resolution is comparable, given the lower resolution setting of 1024 3 .
Apart from CDM, we run cFDM simulations over a range of axion masses m = 10 −22 , 7 × 10 −22 , 2 × 10 −21 eV.DM halos are identified using the friends-of-friends (FoF) algorithm with a standard link-ing length of b = 0.2 × (mean inter-particle separation) (Springel et al. 2001).As a minimum halo resolution, we require all halos to be composed of at least 200 DM resolution elements.We adopt a Planck cosmology (Planck Collaboration et al. 2016) with Ω m = 0.3089, Ω Λ = 0.6911, h = H 0 /100 = 0.6774 and σ 8 = 0.8159.Initial conditions are set up at z = 127, using n s = 0.9665 for the primordial power spectrum of CDM and as input to AxionCamb.We evolve the boxes to redshift z = 3.42.

Estimating HI Fractions
In our simulations, we aim to estimate the individual HI fraction in each gas cell.The IllustrisTNG model uses a modification of the two-phase interstellar medium (ISM) model (Springel & Hernquist 2003, SH03), which assumes star formation occurs only above a certain mass density threshold.Below the threshold, gas physics is determined by hydrodynamics assuming an ideal gas equation of state.Above the threshold, the gas is assumed to have two phases: cold star-forming clouds with a temperature of 1000 K, and hot ionized gas.The SH03 model gives an effective mass-weighted specific energy to a gas cell by averaging over the cold and hot phases, where u eff is the effective specific energy, x is the mass fraction of cold gas in a cell, and u h and u c are the energy per unit mass of the hot and cold gas, respectively.IllustrisTNG uses a modified effective temperature u eff,TNG given by where u 4 is the energy per unit mass corresponding to a temperature of 10 4 K.The HI fraction in each gas cell stored in the snapshots is based on this effective temperature.For star-forming cells, the effective temperature approach underestimates the HI fraction (Villaescusa-Navarro et al. 2018).One can get a better estimate via a direct calculation that does not use the effective temperature.We model the gas similarly as having two phases, and we assume the cold phase is completely comprised of HI, and the hot phase is fully ionized.The cold mass (i.e.HI) fraction x (Springel & Hernquist 2003) can then be expressed analytically as where we define Here t * is the star formation time scale, Λ net is the net cooling rate, ρ is the gas density, u SN is the specific energy corresponding to a supernova temperature of T SN = 5.73 × 10 7 K, and β is the mass fraction in massive stars.Note that we use the IllustrisTNG values of t * = 3.28 Gyr and β = 0.226.Furthermore, we can convert between specific energies and temperatures using the equation where γ is the adiabatic index, k B is the Boltzmann constant, and µ is the mean molecular weight given by Here, m p is the mass of a proton and X H is the hydrogen mass fraction.
A further correction is applied to the HI fraction estimate by subtracting the molecular hydrogen fraction.We adopt the KMT model (McKee & Krumholz 2009), which expresses the mass fraction of molecular hydrogen to neutral hydrogen (both atomic and molecular) f H 2 in star-forming cells as where s is given by and In the above equations, Z is the gas metallicity (in units of solar metallicity), σ d = Z × 10 −21 cm 2 is the cross-section of dust, µ H = 2.3 × 10 −24 g is the mean mass per hydrogen nucleus, and Σ = ρR is the surface density of the gas, where ρ is the gas density and the volume V of the cell can be translated to an effective radius R = (3V/4π) 1/3 .Fig. 1 shows 2D HI mass projections of post-processed ΛCDM and cFDM simulation data.We choose a frequency channel [ν 0 − ∆/2, ν 0 + ∆/2] of width ∆ν = 272 kHz around the central frequency ν 0 = ν 21 /(1 + z), where ν 21 = 1420.406MHz is the rest frequency of the 21 cm line.The channel width corresponds to a comoving width of where r ν is the comoving distance to redshift for the observational frequency ν.We take the same slice of width ∆L along the x-direction from the four hydrodynamic snapshots at z = 4.94, PCS-paint the HI particle data onto a 3D grid (Hand et al. 2018), and project along the x-direction.We see how small-scale structure is visibly suppressed as the axion mass m is reduced from infinity (CDM, top left) down to m = 10 −22 eV (bottom right).

21 cm Cosmology
We focus on the post-reionization era of z < 5.3, in which most of the hydrogen has been ionized and the volume-averaged intergalactic medium (IGM) HI fraction has been reduced to xHI ≈ 10 −4 according to Lyman-α (Lyα) and Lyman-β (Lyβ) effective optical depth measurements (Yang et al. 2020;Bosman et al. 2021).The surviving HI resides in dense regions of matter that are able to selfshield against ionizing background radiation.While most of the HI is found within halos, typically of mass M h ∼ 10 10 − 10 13 M ⊙ , at z = 5 around 12% of the HI is found outside of halos (Villaescusa-Navarro et al. 2018).
The ratio between the number of HI atoms in the ground state and the excited spin state depends on the spin temperature T s , and is given by where n 0 and n 1 are the number of ground and excited spin states respectively and h is Planck's constant.
The differential brightness temperature T b describes the difference between the CMB temperature T γ and the spin temperature T s and is given by where r is a unit vector that is centred on the observer and points in the direction of observation on the sky.The observational frequency ν is related to the redshift z by the standard relation, Eq. (10).We have also defined the optical depth τ 21 across the 21 cm line at redshift z (Pritchard & Loeb 2012;Liu & Shaw 2020), where n H is the hydrogen density, dv || /dr || the gradient of the proper velocity along the line of sight, c the speed of light and A 10 the Einstein spontaneous emission coefficient for the hyperfine transition.
Although observationally the HI can be part of either cold (T ≲ 100 K) or warm (T ≳ 5000 K) neutral medium, both temperatures are warmer than the CMB temperature (which is 19.1 K at z = 6) in the post-reionization era due to heating by sources of radiation such as stars and galaxies.In IllustrisTNG subgrid modeling, lowtemperature cooling channels are not included, leading to a putative cooling floor of 10 4 K.We indeed find that very few gas cells have temperatures below ≈ 3 × 10 3 K in the post-reionization era.Due to Lyα coupling (see e.g.Furlanetto et al. 2006), the spin temperature of HI is thus expected to be higher than the CMB temperature, and we expect to see the 21 cm line in emission.Observationally, the fiducial spin temperature of T s = 100 K is often assumed for DLAs due to a lack of specific measurements, which corresponds to assuming that all of the absorbing gas is in the thermally stable cold neutral medium phase (Murray et al. 2018).In case the gasphase metallicity, dust abundance and background UV field can be estimated, T s can be inferred, but reported values rarely fall below T s = 50 K (Zwaan et al. 2015).It is only in very rare environments where the Lyα radiation field is extremely weak and collisional processes are not effective, that the spin temperature might approach or even fall below the CMB temperature, challenging our assumption T γ /T s ≈ 0.Even though certain DLAs at low redshift z < 1 can have peak optical depths (as inferred from the 21 cm line center) of τ 21 ≈ 0.5 or beyond (Curran 2017;Sadler et al. 2020), it is customary to Taylorexpand the brightness temperature, In fact, the assumption T γ /T s ≈ 0 compels us to Taylor-expand to first order, i.e. these two simplifications that are very common in the literature go hand in hand.We expect the peak values of the brightness temperature field δT b as well as the peak values of the specific intensity field (see Sec. 4) to be slightly overestimated.However, we choose to maintain consistency with the literature since improving the modeling would necessitate a better understanding of the spin temperature T s in optically thick regions.Eq. ( 11) makes clear that the observation of the 21 cm signal can in turn be used as a tracer of matter or as an indirect probe of other properties of our universe, such as its ionization state or temperature.It is also sensitive to cosmology, since the distribution of hydrogen is based on the large-scale distribution of matter, as is the dv || /dr || velocity term, since it includes peculiar velocities in addition to the Hubble flow.Specifically, we make use of the plane-parallel approximation to displace the positions of Voronoi cells from real space (r) to redshift space (s) through where is the familiar Hubble parameter and v || (r) is the peculiar velocity of the cell along the line of sight.While it may seem that it is difficult to probe any of the effects contributing to T b cleanly, different phenomena tend to dominate at different redshifts.Note that in our IllustrisTNG modeling approach, there is a simplifying assumption in the implementation of reionization.The UV background is turned on gradually around z ≈ 10 as a homogeneous, isotropic radiation field, instead of implementing more realistic patchy reionization histories (Byrohl et al. 2021;Molaro et al. 2022;Bird et al. 2023;Puchwein et al. 2023) which significantly af-fect the abundance, distribution, and properties of low column density HI absorbers, changing the shape, distribution, and strength of absorption lines within the Lyα forest.Inhomogeneous reionization not only modifies the post-reionization 21 cm power spectrum and the values for the cosmological parameters inferred from IM experiments (Long et al. 2023) but also the abundance and metallicity distribution of DLAs (Hassan et al. 2020a).We ignore such effects in this work since in the post-reionization era, the assumption of a uniform UV background leads to good agreement with observations of the mean Lyα flux, the HI abundance and the distribution of intermediate and high HI column densities (see Sec. 3 and e.g.Villaescusa-Navarro et al. 2018).
In Fig. 2, we show 21 cm maps in CDM and cFDM cosmologies at z = 4.94 using this approximation.While the brightness temperature δT b is strictly positive throughout the maps, its magnitude spans a large range of over 7 orders of magnitude.The maximum brightness temperature reach δT b,max ≈ 350 mK and are typically associated with DLAs (see Sec. 3.4).Note that δT b,max is highly dependent on numerical resolution.For instance, peak values at z = 4.94 for the 512 3 runs are notably reduced to δT b,max ≈ 95 mK.

HIGlow Modeling of HI
To model HI distributions in the post-reionization era, we utilize normalizing flows, a class of generative models which are hailed for their expressiveness in modeling complex distributions with exact likelihoods and for providing an invertible mapping between a simple base distribution and the target distribution, enabling easy sampling and generation of realistic data.Details of our HIGlow implementation are given in App. A.
The top four rows of Fig. 3 show 64 × 64 HI training maps across CDM and cFDM cosmologies at z = 4.94.The bottom four rows show generated HIGlow HI samples with identical seeds across the cosmologies.The samples visually have very similar features compared to the input training data.The model has learned the effect of changing the axion mass m on the HI maps, and the characteristic suppression of small-scale structure at lower axion mass is apparent.We also show HIGlow samples for a synthetic cosmology corresponding to a new axion mass m = 3 × 10 −22 eV.The generated samples display subtle small-scale features that occupy an intermediate position between those of the m = 10 −22 eV and m = 7 × 10 −22 eV samples, aligning with our expectations.We have conducted comprehensive validation tests, the results of which are outlined in the following section.
To ensure that the statistics of the generated HI maps follow the statistics of the simulation data, we use two validation metrics: HI mass PDFs and HI power spectra.Fig. 4 shows the HI mass PDFs as well as the mean power spectrum and the standard deviation for 1000 simulation data subcube projections (dotted lines) and the same number of HIGlow-generated samples (solid lines).The HI mass PDFs of the generated samples closely follow the target distribution, especially for HI masses below m HI ≈ 10 5 M ⊙ /h.In the range m HI = 10 5 − 10 7 M ⊙ /h, HIGlow captures the distribution less accurately, which is a result of the rarity of such high-mass pixels.For instance, m HI = 10 6 M ⊙ /h pixels are 6 − 7 orders of magnitude less likely than pixels with m HI,peak ≈ 20 M ⊙ /h.The PDF bump around this peak value broadly corresponds to the Lyα forest in the post-reionization era (also see Zamudio-Fernandez et al. 2019).
Note that before training HIGlow, the data undergoes sigmoid normalization as per Eq.(A12).The samples depicted in Fig. 3  presented in sigmoid space, and similarly, as are the power spectra showcased in Fig. 4.This choice enhances the clarity of trends concerning the axion mass m.We find that the target power spectra are reproduced with high fidelity, in particular the mean and standard deviation.The sigmoid-normalized power spectra exhibit the opposite trend with axion mass compared to the 3D non-sigmoidnormalized power spectra (see Sec. 3.3 and e.g.Carucci et al. 2015).We have checked and find that while the projection effect reduces the relative difference between cosmologies as reflected in 2D vs 3D power spectra, the sigmoid normalization reverses the trend.
In the evaluation of HIGlow, a conditional generative model, testing its capability to generate diverse samples spanning the entire latent space of axion masses is crucial.We test this functionality by choosing a new axion mass of m = 3 × 10 −22 eV, unexplored by the simulations, and generating 1000 random images with HIGlow (see Fig. 3).Note that (unlike for CDM and the other three axion masses) we only show the HIGlow generated data for m = 3 × 10 −22 eV, as . Statistical properties of simulation data (dotted lines) and HIGlowgenerated HI samples (solid lines) at z = 4.94 in CDM and cFDM cosmologies.We show HI mass PDFs (top panel), the median of 2D power spectra among 1000 random samples (middle panel) and their standard deviation (bottom panel).In cyan, we add the results for a conditional axion mass of m = 3 × 10 −22 eV (see Fig. 3).This case does not have corresponding simulation data.The Nyquist frequency k Ny = πN 1/3 /L box is shown as a vertical line.
the corresponding simulation data does not exist for this mass.The statistics of the synthetic cosmology are shown in Fig. 4, exhibiting a striking resemblance to the anticipated outcome.Falling in between the curves corresponding to m = 10 −22 eV and m = 7 × 10 −22 eV, the generated distribution effectively showcases the interpolation prowess in latent space.While a direct comparison to actual simu-lations at m = 3 × 10 −22 eV would be needed for a more definitive test, the achieved success can be attributed to the monotonic albeit nonlinear influence of the axion mass m on the distribution of HI.

Overall HI Abundance
Here we investigate the overall HI abundance Ω HI = ρHI (z)/ρ 0 crit relative to the critical density of the Universe today ρ 0 crit , with ρHI (z) being the mean HI density at redshift z.Direct measurements from HIPASS (Barnes et al. 2001) and ALFALFA (Jones et al. 2018;Oman 2022) have been used to detect individual extra-galactic objects, allowing constraints on the HI mass function and abundance of low redshift (z ≈ 0) HI at around Ω HI ∼ (3.9±0.6)×10−4 .At slightly higher redshift, cross-correlations with optical tracers (e.g.DEEP2 or WiggleZ) lead to estimates of the HI abundance at z ∼ 0.8 to be Ω HI b HI ∼ 6.2 +2.3 −1.5 ×10 −4 (Switzer et al. 2013).More indirect measurements of DLAs allow us to constrain the HI abundance up to z ∼ 5. Therefore, in the post-reionization era, we have a well-defined target for the amplitude of the 21 cm signal, which is proportional to the square of the HI abundance, P 21 ∝ Ω 2 HI .We show the value of Ω HI (z) in the four simulated cosmologies in Fig. 5 in the redshift range z = 3.42 − 4.94.HI abundances estimated in redshift space are typically lower than in real space since the 'squashing' effect from coherent large-scale flows is subdominant compared to the Fingers of God effect driven by virial motions (see Sec. 3.3).To compare our predictions to observations, we first note that in a typical DLA survey, HI column densities are estimated from absorption spectra of quasars.The HI frequency distribution or column density distribution function (CDDF, also see Sec. 3.4) is typically written as where N is the number of lines of sight with column densities between N HI and N HI + N HI , and is the so-called absorption distance.The absorption distance has the property that absorbers with non-evolving number density n a and cross-section σ a (both in proper length units) will produce a constant number of absorption lines per unit absorption distance, i.e. dN/dX = const.The distribution f HI (N HI , X) spans ten orders of magnitude from around N HI = 10 12 cm −2 , below which Lyα absorbers produce absorption lines too weak to be recognized, to 10 22 cm −2 above which absorbers are very rare.Quasar absorbers with HI column densities below ≃ 10 17.3 cm −2 are called the Lyα forest, and physically represent the low-density and highly ionized gas that resides in the IGM.Systems with column densities above N HI,DLA = 10 20.3 cm −2 are called DLAs, which typically correspond to extra-galactic regions and are self-shielded against external radiation.Systems with column densities in between these two regimes, are called Lyman Limit Systems (LLSs) in the range 10 17.3 − 10 19.0 cm −2 and sub-DLAs in the range 10 19.0 − 10 20.3 cm −2 .The DLAs account for most of the HI, but it is the Lyα forest that accounts for most of the gas (neutral + ionized).
The total gas density can be inferred from the HI CDDF via where µ = 1.3 is the mean molecular mass of the gas.In the discrete where the sum is calculated for systems with log 10 N HI ≥ 20.3 along lines of sight with a total pathlength of ∆X.As noted in (Zafar et al. 2013;Crighton et al. 2015), we can convert the gas mass in DLAs Ω DLA g to the neutral hydrogen mass via where µ accounts for the mass of helium and δ HI = 1.2 estimates the contribution from systems below the DLA threshold of 10 20.3 cm −2 .In Fig. 5, we add a compilation of measurements, each accompanied by error bars, sourced from multiple studies (Noterdaeme et al. 2009;Songaila & Cowie 2010;Zafar et al. 2013;Crighton et al. 2015;Bird et al. 2017;Ho et al. 2021).Note that all the estimates have been corrected for helium content and cosmological factors.Specifically, both the Hubble constant H 0 and the absorption distance ∆X in Eq. ( 15) have been accounted for.
Within the error bars, the agreement between the results from simulations and observations is good.Hydrodynamic simulations like ours using the IllustrisTNG galaxy formation module typically agree better with observations than semi-analytic models (Lagos et al. 2014) and are comparable to the agreement found using the EAGLE simulations (Rahmati et al. 2015).The trend of decreasing HI abundance towards smaller axion mass m agrees well with the WDM results of Carucci et al. (2015) 1 .
We note that similar to Villaescusa-Navarro (2018), the overall HI abundance depends on resolution.The lower-resolution N = 512 3 simulations contain ∼ 10% less HI in the explored redshift interval z = 3.42 − 4.94 than the higher-resolution N = 1024 3 runs shown in Fig. 5. HI can reside in small halos of mass M h ∼ 10 9 M ⊙ /h, few of which are resolved by the N = 512 3 runs.For halo-based HI estimations such as those based on the halo model (Villaescusa-Navarro et al. 2018), the trend with changing resolution is most intuitive, but it is also reproduced in particle-based models as we find.

HI Power Spectrum and Bias
In Fig. 6, we show the total matter and HI power spectrum estimated using a third-order piecewise cubic spline (PCS) mass-assignment scheme (MAS, Hand et al. 2018) and subsequent Fourier transformation, with X = {tot, HI} and N modes being the number of modes lying in the spherical shell k of width δk = 2π/L box .When calculating the HI power spectrum, we do not account for the fact that the actual density profile of the gas particles is not given by a uniform cube but described by the SPH kernel.A reliable amplitude mode correction implementation is beyond the scope of this paper (also see Villaescusa-Navarro et al. 2014).However, we compensate for the MAS window function to improve the power spectrum fidelity on scales close to the Nyquist frequency k Ny = πN 1/3 /L box , where N 1/3 is the number of cells per box side.
In cFDM cosmologies, the total matter power spectrum in Fig. 6 follows the well-known trend of small-scale suppression for wavenumbers above k 1/2 (Marsh 2016).The effect of peculiar velocities, also called redshift-space distortions (RSDs), can be clearly discerned.On large scales, the clustering of matter in redshift space is enhanced due to the Kaiser effect (Kaiser 1987).On small scales, the peculiar velocities, particularly inside halos, give rise to the Fingers of God, suppressing the amplitude of the total matter power spectrum (Tegmark et al. 2004).
HI is more clustered than the total matter since the UV background ionizes hydrogen in environments that are not dense enough to self-shield.As opposed to the suppression beyond k 1/2 in the total matter clustering, the HI power spectrum does not exhibit a smallscale cut-off since most of the HI is trapped in intermediate-and high-mass halos.Contrary to a naive expectation, HI clustering increases with decreasing axion mass m.Since most of the HI is locked inside DM halos in the post-reionization era, the higher halo bias in cFDM cosmologies compared to CDM (see e.g.Dunstan et al. 2011) leads to an increased HI clustering.This is in agreement with Carucci et al. (2015) who showed that the amount of HI per total mass in a given DM halo mass bin is very similar between CDM and WDM cosmologies for particle-based HI modeling approaches.As a result, the spatial distribution of HI is more biased in the cFDM models than in CDM, with the bias increasing with decreasing axion mass m.In the extreme cFDM model with m = 10 −22 eV, deviations from CDM in the HI power spectrum can reach more than 300% at z = 4.94 (in agreement with WDM results from Carucci et al. 2015).
The HI bias b HI (k) = P HI (k) describes the bias between the distributions of HI and the to- HI Bias tal matter.While the HI bias defined using Eq. ( 19) suffers from stochasticity (shot noise), it is closer to observations than results inferred using the alternative definition involving the HI-matter crosscorrelation, b HI (k) = P HI-tot (k)/P tot (k) (Sarkar et al. 2016;Castorina & Villaescusa-Navarro 2017).As shown in Fig. 6, the large-scale halo bias in CDM at z = 4.94 attains values b HI ≈ 2.5 in redshift space and b HI ≈ 3.5 in real space, in agreement with previous work (Sarkar et al. 2016;Villaescusa-Navarro et al. 2018;Wang et al. 2021).On the smallest (strongly non-linear) scales probed, it increases up to b HI ≈ 20.In cFDM cosmologies, the HI bias increases with decreasing axion mass m, and for the extreme m = 10 −22 eV model, the bias ratio between cFDM and CDM can reach values above ≈ 2.5.
The clustering of the 21 cm signal which is modeled following Sec.2.3 is shown in Fig. 7. Its magnitude ∝ Ω 2 HI (z) in the postreionization era is several orders of magnitude lower than before and during reionization, and exhibits a clear trend with decreasing axion mass m.The SKA1-Low system noise estimated using 21cmSense (Pober et al. 2013(Pober et al. , 2014) ) is shown for instrument parameters from Table 3, assuming an optimistic foreground wedge model in which all k modes inside the primary field of view are excluded (Pober et al. 2014).Note that we exclude the error arising from sample variance, which predominantly affects larger scales.We find that at z = 4.94, the 21 cm power spectrum will be detected by this telescope up to scales k ≈ 9 h/Mpc in CDM and the less extreme cFDM models m ∼ 10 −21 eV, while it can be detected up to k ≈ 15 h/Mpc in the more extreme m = 10 −22 eV model.
Since we assume that T γ /T s ≈ 0 in our modeling, the relative difference in the 21 cm power spectrum between cFDM and CDM is identical to the relative difference in the HI power spectrum shown in Fig. 6.We also show the normalized error ∆ 2 noise /∆ 2 21,CDM for SKA1-Low, representing the 1σ system noise error on the 21 cm dimensionless power spectrum of the model with CDM.While the more extreme m = 10 −22 eV model can be distinguished from CDM at the 2σ confidence level across all scales (k < 80 h/Mpc) resolved by the simulations, SKA1-Low will not be able to discriminate among CDM and the less extreme cFDM models m ∼ 10 −21 eV with 1080 hours of observations on small scales k ≈ 70 h/Mpc.

HI Column Density Distributions
Another quantity commonly employed to study the distribution of HI in the post-reionization era is the HI CDDF, Eq. ( 13).As mentioned in Sec.2.2, HI column densities of Lyα clouds, LSS, sub-DLAs and DLAs are inferred observationally from quasar spectra.
To compare mock HI CDDFs to observational datasets, we first assign the HI to gas particles following Sec.2.2.Then, the value of the column density N HI along an arbitrary line of sight (LOS) can be computed following a method from Villaescusa-Navarro et al. ( 2014) which we briefly summarize here.We first project all the Figure 7. 21 cm power spectrum in CDM and cFDM cosmologies at z = 4.94.We show results for the 'dimensionless' power spectrum ∆ 2 21 (k) = k 3 P 21 (k)/(2π 2 ), in both real space (dotted lines) and redshift space (solid lines).Black crosses denote the expected SKA1-Low 1σ system noise ∆ 2 noise estimated using 21cmSense (Pober et al. 2013(Pober et al. , 2014) ) for 1080 hours of observation in a 8 MHz bandwidth with 290 frequency channels as per Table 3, for an optimistic foreground wedge model.The vertical dashed line denotes the Nyquist frequency k Ny = πN 1/3 /L box (in black).In the bottom panel, we present the relative difference between cFDM and CDM, identical to the relative difference in the HI power spectrum, see Fig. 6.The normalized error ∆ 2 noise /∆ 2 21,CDM on the CDM 21 cm power spectrum for the same instrument parameters in SKA1-Low is shown by a shaded region.21 cm clustering increases for decreasing axion mass m.
gas particle positions onto the XY Cartesian plane, draw LOSs on a regular grid perpendicular to the XY plane from Z = 0 to Z = L box .For a given LOS, if the minimum distance between any gas particle i and the LOS (i.e.impact parameter b i ) is smaller than the smoothing length h i of the gas particle, then we pick up a contribution based on the integrated density along the path that the LOS intersects the physical size of the gas particle: where N HI,i is the column density due to particle i, having HI mass m HI and SPH smoothing length h i while m H is the mass of the hydrogen atom.The relation between the radius r and the integration variable l is given by r 2 = b 2 + l 2 , with l 2 max = h 2 i − b 2 .Our regular grid consists of 10, 000 × 10, 000 points, i.e. the number of LOSs  2021) perform a DLA analysis using SDSS DR12 and DR16Q data, respectively.The redshift ranges in the legend refer to the (sub-)DLA redshifts and not the emission redshifts of the quasars. is 10 8 .In view of our small box size L box = 40 h −1 Mpc, the probability of encountering more than a single absorber with a large column density ∼ 10 19 cm −2 along the LOS is negligible.Our results are thus converged for sub-DLAs and DLAs for which we can, as a result, safely ignore light-cone effects.We repeated the tests mentioned in Villaescusa-Navarro et al. (2014) to verify that the grid is fine enough to achieve convergence.For instance, when slicing the box into slabs of different widths and computing the column densities from any of those slabs, the frequency distribution does not change above ∼ 10 19 cm −2 .
The inferred frequency distribution of HI column densities in CDM and cFDM cosmologies is shown in Fig. 8.We assume a weak redshift dependence of the HI column distribution in the range z = 2 − 4, which can be understood in the modeling framework of Theuns (2021) (Ho et al. 2021), the trend beyond z > 4 is in fact still unclear due to statistical uncertainties in the measurements.To better align with the observations, the results in the CDM and cFDM models are shown at z = 3.42 rather than z = 4.94 as for the bulk of this paper.
All simulated cosmologies show good agreement with the observations in the sub-DLA and DLA regimes.The strongest DLA absorbers above N HI = 3 × 10 22 cm −2 are not modeled well in any Table 1.Best-fit parameters M 0 in Eq. ( 22) for CDM, m = 2 × 10 −21 eV and m = 7 × 10 −22 eV cosmologies and M break , β in Eq. ( 23) for the m = 10 −22 eV cosmology across redshifts z = 3.42 − 4.94. of the simulated cosmologies with the predicted CDDF undercutting the measurements, which can be traced to our moderate box size of L box = 40 h −1 Mpc.The cFDM cosmologies start to diverge from the CDM result in the sub-DLA regime.As a result, with larger data sets and better mean flux measurements, one might be able to put competitive astrophysical constraints on FDM models using sub-DLA surveys.In the LLS and the Lyα forest regimes, the abundance of absorbers is expected to differ significantly between CDM and cFDM cosmologies.However, in these regimes we can neither estimate the frequency distribution reliably (as mentioned above) nor is the assumption of a homogeneous UV background valid.

DLA Cross-Sections
Next, we focus on absorbers with high column densities above N HI,DLA = 10 20.3 cm −2 , i.e.DLAs and their cross-sections.The latter are of interest since the observed rate of incidence of DLAs tells us the product of the number density of DLAs times their cross-section (Font-Ribera et al. 2012;Pérez-Ràfols et al. 2018).Both the number density and the cross-section of absorbers is expected to be sensitive to DM physics.In cFDM cosmologies, the high column density of cosmic filaments at z = 2 − 5 (Gao et al. 2015) should be reflected in the cross-section distribution, which we now investigate for the first time.
We compute the DLA cross-section as follows.For each DM halo, we select the gas particles belonging to the halo and we throw random LOSs within the halo virial radius.After computing the column density for each LOS, the DLA cross-section is readily obtained as where R vir is the halo virial radius, n DLA is the number of LOSs with column density above 10 20.3 cm −2 and n tot is the total number of LOSs used.The virial radius is determined as per the Bryan & Norman (1998) spherical collapse result, but we checked that our results do not change if we use a different definition for R vir .We use 40, 000 random LOSs for each DM halo, and find that our results do not change significantly when we increase the number.We show results in Fig. 9.We find that the mean DLA crosssection increases with halo mass, which apart from the m = 10 −22 eV cosmology is well fitted by the following function In agreement with the CDM analysis of Villaescusa-Navarro et al.
(2018), we find that the slope of the cross-section for large halo masses has best-fit value α = 0.82 and the characteristic halo mass where the DLA cross-section decreases exponentially has a strong correlation with the overall normalization of the function, The cutoff exponent is well approximated by β = 0.85 • log 10 (N HI,DLA /cm −2 ) − 16.35, the only redshift dependence remaining in M 0 and A. At z = 4.94 which is shown in Fig. 9, the characteristic halo mass is M 0 = 1.6 × 10 9 h −1 M ⊙ in all three cosmologies.In Table 1, we present best-fit values for M 0 across redshifts z = 3.42 − 4.94, which decreases with redshift implying that less massive halos can host HI at higher redshift (also see Villaescusa-Navarro et al. 2018).While a small secondary peak emerges in the cross-section distribution for the m = 7 × 10 −22 eV model, the m = 10 −22 eV cosmology exhibits a bimodal distribution at z = 4.94.The median crosssection is significantly higher than in the other three cosmologies for halo masses below M ≈ 3 × 10 9 M ⊙ /h, which is about an order of magnitude below the half-mode mass M 1/2 (Marsh 2016).We thus confirm the prediction of Gao et al. (2015) that the high column density of cosmic filaments in cFDM/WDM cosmologies is reflected in the properties of DLAs.This is also in agreement with cosmic web studies at high redshift, which show that the mean DM overdensity in filaments at z ≈ 5 is around 30% higher in the m = 10 −22 eV cosmology than for CDM (Dome et al. 2023b).The power law + cutoff parametrization of Eq. ( 22) thus gives a poor fit for the m = 10 −22 eV model, and instead we adopt a double power law The dependence of σ DLA on cosmology at the high-mass end (i.e.M > M break ) is so weak that the normalization A can still be well approximated by A•M 0 = 14100 h −3 kpc 2 M ⊙ , where M 0 is taken from Eq. ( 22).The normalization at the low-mass end B is determined by ensuring that the two power laws intersect at M break .We provide the best-fit values for M break and β in Table 1.
Using the mean DLA cross-section, we can estimate the DLA bias as where n(M, z) denotes the halo mass function (HMF) and b(M, z) the halo bias.We take the HMF from Sheth & Tormen (2002) for CDM and apply the Schneider et al. (2012) rescaling for cFDM.For the halo bias, we use the fitting formula in Tinker et al. (2010), which expresses the bias b(M, z) as a function of the peak height ν = δ c /σ(M), where δ c = 1.686 is the spherical collapse threshold and σ 2 (M) the linear matter variance.This parametrization provides a reliable fit for the bias of both CDM and cFDM halos when adopting the variance σ 2 (M) in the respective cosmology (Dunstan et al. 2011).
In Table 2, we show the estimated values for the DLA bias in the various cosmologies across redshifts z = 3.42 − 4.94, for absorbers with column density N HI > N HI,DLA = 10 20.3 cm −2 .We find that the DLA bias increases with redshift at the explored post-reionization redshifts.Since the gas in the IGM is denser and the amplitude of the UV background decreases with redshift, this trend is expected for the DLA bias as well as the HI bias.Typically, both biases have similar magnitudes, b DLA ≃ b HI , following the theoretical arguments in Castorina & Villaescusa-Navarro (2017) and comparing the values in Table 2 to the b HI estimates from Fig. 6 at z = 4.94.In cFDM cosmologies, we find the redshift-independent trend that the DLA bias increases with decreasing axion mass m, from b DLA = 2.3 for CDM to b DLA = 3.8 for the m = 10 −22 eV model at z = 4.94.This trend is primarily a result of the increased halo bias in cFDM cosmologies, as discussed in Sec.3.3.
It is useful to compare our estimates of the DLA bias to observations.Through cross-correlating the DLAs with the Lyα forest, Pérez-Ràfols et al. ( 2018) measured a DLA bias factor of b DLA = 1.99 ± 0.11, while using cross-correlations with CMB lensing data Alonso et al. (2018) and Lin et al. (2020) measured b DLA = 2.6 ± 0.9 and b DLA = 1.37 +1.30  −0.92 , respectively.All three estimates are for a median DLA redshift of z median = 2.3 and agree well with our CDM estimates from Table 2, if we assume a linear relationship from z = 3.9 to z = 2.3 (see e.g.Villaescusa-Navarro et al. 2018).With more powerful (combined) datasets and a better understanding of the theoretical systematics involved in our use of Eq. ( 24), the observed DLA bias could be used to put constraints on FDM-like cosmologies.a For imaging: N chan = 29; for 21 cm power spectra: N chan = 290.
Table 3. Summary of SKA1-Low instrument parameters used in this study to obtain the thermal noise sensitivity.The FWHM of the primary beam agrees well with Kakiichi et al. (2017) and Hassan et al. (2020b).Natural sensitivities A eff /T sys as a function of frequency are tabulated in Braun et al. (2019).

MOCK RADIO MAPS FOR SKA-LOW
Here we investigate the prospects of imaging the HI distribution in the post-reionization era using SKA1-Low.Since SKA maps and summary statistics such as the 21 cm power spectrum are poised to shed light on the nature of DM (Carucci et al. 2015;Bauer et al. 2021), we study how ML models such as normalizing flows may support such efforts.While in an average region of the Universe, the 21 cm signal may not be strong enough to resolve in imaging, regions with strong HI concentrations are more likely to be resolved.We thus focus on whether SKA1-Low might be able to detect a few individual bright HI peaks and how these peaks link to DM physics.We begin by creating brightness temperature maps from our simulated distribution of HI following Secs.2.2 and 2.3.We adopt a channel width of ∆ν = 272 kHz, corresponding to a subbox with side length 2.5 h −1 Mpc in agreement with the HIGlow modeling of Sec.3.1.We relate the brightness temperature excess to the specific intensity using the relation and project the I ν (r, ν) field onto a regular grid of resolution 1024 × 1024 in the XY plane.

Angular Resolution
The intrinsic resolution of our simulation ∆x = L box /N ≈ 40 h −1 kpc corresponds to angular scales where D c (z) is the comoving distance to redshift z = 4.94.On the other hand, the angular resolution of the radio interferometric observation is determined by the maximum baseline B max via θ A = λ/B max , where λ is the signal wavelength at the observation redshift z.Since the maximum baseline of SKA1-Low is B max = 65.4 km corresponding to ≈ 4.0 ′′ , SKA1-Low could get close to such We account for the chosen instrumental resolution by smoothing the I ν field with a Gaussian window of angular radius θ A , with θ A being the FWHM of the synthesized beam given in Table 3. Radio interferometers are not sensitive to the mean value of the specific intensity and can only measure deviations from the mean, hence we recenter the specific intensity field by making the transformation I ν → I ν − ⟨I ν ⟩.Finally, we multiply the I ν map with the beam solid angle to arrive at the total flux within a single beam, i.e. the resulting image maps will be in units of flux per beam.We repeat the procedure until we find the slice from the simulation box that contains the highest value of I ν , and report results in Table 4.

Thermal Noise
To generate noise maps and to estimate the noise levels we largely follow Hassan et al. (2019).Every pair of stations in SKA1-Low will record noise along with the visibilities in the uv-plane.This noise is uncorrelated between measurements and can be represented by Gaussian random numbers of mean zero and standard deviation (e.g.Ghara et al. 2017) where t int is the integration time to observe a single visibility at frequency resolution ∆ν, A eff is the effective area of each station and T sys is the system temperature which combines the receiver temperature and the sky temperature.Their combination A eff /T sys , the natural sensitivity, is frequency dependent and we retrieve the most up-to-date tables from Braun et al. (2019).Eq. ( 28) provides the flux error for a single baseline which we convert to Kelvin using where Ω = π/(4 ln 2)θ 2 max is the full beam area.We present rms noise errors for a single baseline in Table 3, along with the SKA1-Low instrument parameters.
To generate a 2D thermal noise realization we follow four steps: (i) We first generate Gaussian random noise (both the real and imaginary parts) with zero mean and rms ⟨|N| 2 ⟩ from Eq. ( 28) in the N pix × N pix grid in Fourier space, where N pix = θ max /θ A .(ii) We compute the daily uv coverage N uv which represents the total number of baselines that observe a given uv pixel using 21cmSense .
We account for the presence of multiple baselines in a uv grid point by scaling the noise in the (i, j)th pixel by a factor 1/ N i j uv at nonzero uv coverage pixels.(iii) By averaging over the total number of days observed, we reduce the noise further by a factor of 1/ N days .(iv) We obtain the real space noise map by inverse Fourier transforming the noise from Fourier space at each frequency channel.
The rms noise values for different values of the angular resolution θ A obtained using the above four-step procedure are given in Table 4.To determine whether there will be enough sensitivity in the maps to identify the cosmological HI using SKA1-Low, we calculate the SNR values of the peak fluxes.While the SNR for θ A = 2 ′ is modest (SNR > 3), for θ A = 1 ′ the SNR is very low (SNR < 2).We thus conclude that the imaging prospects of SKA1-Low for 1080 observing hours and other fiducial parameters listed in Table 3 are moderate at θ A = 2 ′ , with a rapidly declining SNR for lower values of θ A .Note that some of the previous works assume a higher number of antenna stations in the SKA1-Low design (Villaescusa-Navarro et al. 2014;Ghara et al. 2017;Hassan et al. 2019Hassan et al. , 2020b)), corresponding to a higher number of baselines and thus lower thermal noise.In particular, for that reason Villaescusa-Navarro et al. ( 2014) arrived at higher peak SNR estimates for SKA1-Low at z = 4. Accounting for angular resolution and thermal noise, in Fig. 10 we show mock radio maps at z = 4.94 in CDM for two different angular resolutions, θ A = 1 ′ and θ A = 2 ′ .Notably, a prominent peak (SNR= 3.25) can be observed in the lower half for θ A = 2 ′ , corresponding to fluxes per beam around 2µJy.However, enhancing the angular resolution comes at a trade-off of increased thermal noise.This effect is depicted in the top panel of Fig. 10.Here, when adopting θ A = 1 ′ , the peak becomes more challenging to distinguish (SNR=1.81).

HIGlow Interpolation
The imaging prospects for cFDM cosmologies are similar to CDM, see peak flux and SNR values in Table 4.At our fiducial redshift of z = 4.94, we find a trend of decreasing peak flux I max for lower values of the axion mass m.We could not verify this trend at all redshifts and thus do not believe that it is of physical origin.However, as a proof of concept, we take advantage of this trend and demon- strate the ability of HIGlow to interpolate the latent space of axion masses.
As mentioned in Hassan et al. (2022), assuming that HIGlow correctly captures the conditional distribution, it can be used to generate new samples for conditional parameters on which the model has not explicitly been trained.To this end, we generate 1000 HIGlow samples for a synthetic cosmology with m = 3 × 10 −22 eV (see Sec. 3.1) and quote the maximum peak flux I max among all samples in Table 4.We find that the values are in between the m = 7 × 10 −22 eV and m = 10 −22 eV values inferred from the full-box approach, indicating the success of HIGlow to interpolate its conditional latent space after training.

CONCLUSIONS
21 cm intensity experiments with future radio telescopes like the SKA offer a promising source of new astrophysical and cosmological information about the post-reionization era.As such, we can expect competitive constraints on the DM particle mass m and noncold DM fraction (see e.g.Giri & Schneider 2022).In this work, we model the distribution of neutral hydrogen (HI) using highresolution hydrodynamical N-body simulations in CDM and three instances of cFDM cosmologies with axion masses m = 10 −22 , 7 × 10 −22 , 2 × 10 −21 eV, adopting the IllustrisTNG galaxy formation module.Our focus is on the redshift range z = 3.42 − 4.94.We demonstrate the ability of the emerging ML technique of normalizing flows to learn the relevant features in the HI distribution and illustrate the detectability of the brightest HI peaks using SKA1-Low.We summarize our main conclusions as follows: (a) We observe a trend of decreasing HI abundance Ω HI for smaller axion mass m.Comparing to a range of recent measurements, we show that both CDM and the less extreme cFDM models m ∼ 10 −21 eV are in good agreement therewith.However, the statistical and systematic errors on the measurements are too large to put competitive constraints on FDM using Ω HI directly.(b) The HI bias b 2 HI (k) = P HI (k)/P tot (k) as inferred from the ratio of the HI and total matter power spectra increases with decreasing axion mass, as previously establish.We identify that it is the higher halo bias in cFDM cosmologies (e.g.Dunstan et al. 2011) that leads to increased HI clustering.This is also reflected in the amplitude of the 21 cm power spectrum which is up to 300% higher for the m = 10 −22 eV model than for CDM at z = 4.94, similar to (Carucci et al. 2015).We find that SKA1-Low will be able to discriminate among the m = 10 −22 eV model and CDM at the 2σ confidence level on all scales k < 80 h/Mpc resolved by the simulations.(c) The weak dependence of the HI frequency distribution f HI (N HI , X) on redshift in the post-reionization era allows to compare to a range of DLA and sub-DLA measurements covering z = [1.5 − 5.0].These show good agreement with predictions from all our investigated CDM and cFDM models.We show that in the sub-DLA regime, cFDM cosmologies start to exhibit a lower abundance of HI absorbers compared to CDM.(d) DLA cross-sections exhibit a simple power law + low-mass exponential cutoff trend for cFDM models m ∼ 10 −21 eV that is well known for CDM models (Villaescusa-Navarro et al. 2018).For more extreme models with m ∼ 10 −22 eV, we show that a double power law parametrization provides a better fit.The high median crosssection at the low-mass end can be traced to the high column density of cosmic filaments (Gao et al. 2015).Using the mean DLA cross-section, we estimate the DLA bias b DLA and recover the wellknown result of increasing DLA bias with increasing redshift.In addition, we find the redshift-independent trend that the DLA bias  2010), to capture intricate structures resulting from non-linear physics.We adapt the HIGlow framework (Hassan et al. 2022) to facilitate the conditional generation of HI maps with varying axion masses.We use two validation metrics, HI mass PDFs and HI power spectra.As proof of concept, we showcase the ability of HIGlow to interpolate its latent space of axion masses to predict the peak flux value for a synthetic cosmology with m = 3 × 10 −22 eV.
As we advance into the systematics-dominated era and witness the construction of PB-scale radio telescopes such as SKA, in the quest for DM it is imperative to improve the modeling of HI distributions in various CDM and non-CDM cosmologies.In the postreionization era that we focus on, it is easier to disentangle DM signatures from astrophysical processes, which are poorly understood during the epoch of reionization and cosmic dawn.Observations during (e.g.Castellano et al. 2023) and after reionization (e.g.Iršič et al. 2023) are thus complementary and will increase the robustness of DM constraints.However, a more comprehensive analysis of astrophysical processes is needed even in the post reionization era to swiftly explore a wide array of galactic and stellar feedback effects that can be marginalized over.

DATA AVAILABILITY
The IllustrisTNG galaxy formation module specifications are publicly accessible at https://www.tng-project.org/.Post-processing scripts are made available upon reasonable request.
Table A1.Three layers making up a conditional flow block.g s , g t are simple neural networks, and f s , f t are some simple well behaved functions.In particular, a sigmoid function is used for f s because a requirement of the affine coupling block is that this function should not take a value of zero (Behrmann et al. 2020).Note how the model is made conditional by introducing the conditional parameter c into the affine injector and affine coupling layers.
measures a statistical distance between distributions.For a large number of samples from the true distribution, we get the following estimate for the loss function (Papamakarios et al. 2021): In Eq. (A7), the first term is the standard log-likelihood of a multivariate Gaussian.Since we choose simple functions to make up the flow, the Jacobians and their determinants in the second term have simple analytical expressions.Now another great advantage of normalizing flows becomes evident: the exact likelihoods are known and worked with throughout the training.By contrast, other generative models like generative adversarial networks (GANs) and variational autoencoders (VAEs) do not allow us to extract the exact likelihoods.By inverting the flow, normalizing flows thus provide a tractable high-dimensional likelihood for constraint forecasting and parameter inference in astrophysics and cosmology (e.g.Zamudio-Fernandez et al. 2019;Bevins et al. 2022;Hassan et al. 2022).
Our model is naturally a generative one.Once trained, we can start with a Gaussian noise value y, the axion mass parameter c and generate a sample HI map x by applying the inverse flow g: x = g(y, c) = f −1 (y, c). (A8)

A2 Architecture
The challenge of designing a normalizing flow architecture is choosing the functions f i that make up the flow.There are various different existing architectures such as RealNVP (Dinh et al. 2016) and NICE (Dinh et al. 2014).However, we use a Glow architecture (Kingma & Dhariwal 2018) which has been modified to be conditional (Friedman & Hassan 2022).The building block of the unconditional Glow model is a flow step, which is the composition of an actnorm layer, an invertible 1 × 1 convolution, and an affine coupling layer.Multiple flow steps are composed to make a flow block, and multiple flow blocks make up the whole model.The flow blocks are connected using a multi-scale architecture (Dinh et al. 2016) which splits the data in between flow blocks.This approach allows the use of the spatial structure of the HI maps during training.To make Glow conditional, the actnorm and affine coupling layers are modified in a simple way

Figure 1 .
Figure 1.HI mass maps in CDM and cFDM cosmologies at the post-reionization redshift of z = 4.94.The width of the projected slice ∆L ≈ 2.5 h −1 Mpc corresponds to a channel width of ∆ν = 272 kHz.Redshift-space distortions (RSDs) are included by displacing Voronoi cells following Eq.(12).

Figure 2 .
Figure 2. 21 cm maps in various cosmologies at z = 4.94.The comoving projections are over the same width as in Fig. 1.The highest 21 cm brightness temperatures δT b ≈ 350 mK are found in DLAs which can self-shield against external UV radiation.
Figure 3. Simulated data training samples (top four rows) and generated samples with HIGlow (bottom four rows) in the various cosmologies at z = 4.94.The 64 2 maps are 2.5 h −1 Mpc a side.We add HIGlow samples for a conditional axion mass m = 3 × 10 −22 eV.This case does not have corresponding simulation data, illustrated by crosses.The axion mass m decreases from right to left, with increased small-scale suppression.

Figure 5 .
Figure 5.Comparison of mock and observed cosmic HI density estimates.The abundance parameter Ω HI = ρHI (z)/ρ 0 c from simulations in the redshift range z = 3.42 − 4.94 for various cosmologies is shown by dotted lines (in real space) and solid lines (in redshift space), obtained using the procedure outlined in Sec.2.2.Observational measurements are displayed with 1σ error bars: Songaila & Cowie (2010) provide DLA measurements from Keck data; Zafar et al. (2013) quote combined DLA and sub-DLA measurements from ESO UVES; Crighton et al. (2015) quote results from a Gemini GMOS study of DLAs; Noterdaeme et al. (2009), Bird et al. (2017) and Ho et al. (2021) are DLA analyses using SDSS DR7, DR12 and DR16Q, respectively.

Figure 6 .
Figure6.Total matter (left panel), HI power spectrum (middle panel) and HI bias (right panel) in CDM and cFDM cosmologies at z = 4.94.Each panel shows results in real space (dotted lines) and redshift space (solid lines).The vertical dashed lines mark the cFDM half-mode scale k 1/2 (colored) as perMarsh (2016) and the Nyquist frequency k Ny = πN 1/3 /L box (in black).In the bottom row, we present the ratio of results between cFDM and CDM, with the first and third panels displaying the absolute ratio and the second panel showing the relative difference in percent.Both HI clustering and HI bias b 2 HI (k) = P HI (k)/P tot (k) increase for decreasing axion mass m.

Figure 8 .
Figure 8. Column density distribution function of HI absorbers.Results at z = 3.42 in CDM and cFDM cosmologies are traced by solid lines.Data from observations are shown with 1σ error bars: Péroux et al. (2005) quote (sub-)DLA measurements from ESO UVES; Noterdaeme et al. (2012) perform a DLA analysis using SDSS DR9 data; Zafar et al. (2013) quote sub-DLA measurements from ESO UVES; Crighton et al. (2015) quote results from a Gemini GMOS study of DLAs; Bird et al. (2017) and Ho et al. (2021) perform a DLA analysis using SDSS DR12 and DR16Q data, respectively.The redshift ranges in the legend refer to the (sub-)DLA redshifts and not the emission redshifts of the quasars.
based on cosmological accretion of gas onto DM halos.Hence, we can overlay results from observations at redshifts z = [1.5 − 4.0].Since the (sub-)DLAs at z > 4 are rare and typically inferred from noisy spectra, they mostly contribute to enlarged error bars and we can extend the redshift range to z = [1.5 − 5.0].We thus add measurements from Péroux et al. (2005); Noterdaeme et al. (2012); Zafar et al. (2013); Crighton et al. (2015); Bird et al. (2017); Ho et al. (2021) with errorbars to Fig. 8.As shown in the most recent measurement paper

Figure 9 .
Figure 9. DLA cross-section of halos at z = 4.94.The color of each hexbin is proportional to the number of halos.For the CDM (top left), m = 2 × 10 −21 eV (top right) and m = 7 × 10 −22 eV (bottom left) cosmologies, we fit a power law + exponential function (22) while for the m = 10 −22 eV (bottom right) cosmology we fit a double power law, Eq. (23).The vertical dashed line denotes the half-mode mass M 1/2 (Marsh 2016).The small secondary peak in the m = 7 × 10 −22 eV model morphs into a bimodal distribution in the m = 10 −22 eV cosmology.

Figure 10 .
Figure10.Mock radio maps for CDM with a frequency channel of width ∆ν = 272 kHz and synthesized beam FWHM θ A = 1 ′ (top row) and θ A = 2 ′ (bottom row).Left column: Specific intensity field I ν (r) for N pix = 1024 corresponding to our simulation resolution.Right column: Specific intensity field overlaid with instrumental noise for N pix = 12 and N pix = 25, respectively.We construct the noise map following a four-step procedure outlined in Sec.4.2.
increases with decreasing axion mass m, from b DLA = 2.3 for CDM to b DLA = 3.8 for the m = 10 −22 eV model at z = 4.94.This trend is in agreement with Point (b) and is again a result of the increased halo bias in cFDM cosmologies.(e) The prospects of imaging the brightest HI peaks with SKA1-Low at the fiducial redshift of z = 4.94 are moderate (SNR > 3) for angular resolutions θ A ≈ 2, with a rapidly declining SNR for lower values of θ A .Our peak SNR estimates are lower than for the comparable study of (Villaescusa-Navarro et al. 2014) who assumed 911 stations for SKA1-Low as opposed to 512.(f) We demonstrate the ability of normalizing flows, a generative ML model introduced by Agnelli et al. (

)
Hence our loss function is just the negative log-likelihood.Applying the standard change of variables formula for PDFs, we get log p θ (x|c) = log p θ (y) + log det dy d(x, c) = log p θ (y)

Table 4 .
Prospects of detecting the brightest HI peaks using SKA1-Low at z = 4.94.We provide the peak flux I max , noise rms ⟨|N| 2 ⟩ (both in µJy/beam) and SNR values in CDM and cFDM cosmologies for channel width ∆ν = 272 kHz, corresponding to a subbox with side length 2.5 h −1 Mpc in agreement with the HIGlow modeling of Sec.3.1.The peak flux values for the conditional mass of m = 3 × 10 −22 eV (not part of our simulation suite) are obtained by generating 1000 random samples using HIGlow and searching for the strongest peak among all samples.However, due to the compact core layout most of the sensitivity is at the shorter baselines, hence the noise level is lower if one reduces the resolution of the images.In this work, we choose angular resolutions θ A = 0.5 ′ , 1 ′ , 2 ′ corresponding to baselines B max = 8.6, 4.3, 2.2 km.
(Giri & Schneider 2022et al. 2021ons realizing thousands of astrophysical models(Villaescusa- Navarro et al. 2021) and seminumerical modeling techniques(Giri & Schneider 2022) signify a major stride in the pursuit of this goal.It is a pleasure to thank Oscar O'Hara for enriching conversations on the intricacies of radio interferometers, particularly in the context of the future Square Kilometer Array.We express our sincere gratitude to Roy Friedman for his support in bolstering the numerical stability of the HIGlow model.We also appreciate insightful discussions with Emily Wickens.TD acknowledges support from the Isaac Newton Studentship and the Science and Technology Facilities Council (STFC) under grant number ST/V50659X/1.AF is supported by the Royal Society University Research Fellowship.The simulations were performed under DiRAC project number ACSP253 using the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility.The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1.DiRAC is part of the National e-Infrastructure.