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

Tibor Dome1,2, Rumail Azhar1, Anastasia Fialkov1,2

1Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge,
CB3 0HA, UK
2Kavli Institute for Cosmology, Madingley Road, Cambridge, CB3 0HA, UK E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

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 and the HI column
density distribution function decrease 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
Ly${\alpha}$ (DLA) bias which we infer from the cross-section of DLAs. The DLA
cross-section distribution 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.
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.

###### keywords:

methods: numerical - galaxies: high-redshift, intergalactic medium -
cosmology: dark matter, large-scale structure of Universe - radio continuum:
general

††pubyear: 2023††pagerange: Modeling Post-Reionization HI Distributions in
Fuzzy Dark Matter Cosmologies Using Conditional Normalizing Flows–References

##  1 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 observed frequency-redshift
relation $\nu=\nu_{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 $\Lambda$CDM model - faces some challenges
including the lack of direct detection of DM particles and small-scale
observations that are at odds with $\Lambda$CDM predictions (Weinberg et al.,
2015; Bullock & Boylan-Kolchin, 2017). While attempts at unified baryonic
solutions to the $\Lambda$CDM small-scale 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 $\sim
10^{-22}$ eV called axions, corresponding to de Broglie wavelengths of
$\lambda_{\text{dB}}\sim 1$ kpc. The motivations behind considering multiple
species of light axions is vast, encompassing well-established predictions of
string/M-theory (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 $\lambda_{\text{db,c}}$
which scales as $\lambda_{\text{db,c}}\sim(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
$\Lambda$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\times 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., 2022). 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. (2018, 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 slightly modified version (see Sec. 2.4) of the
conditional normalizing flow framework HIGlow (Friedman & Hassan, 2022) and
show that the efficiency of the model in sample generation, coupled with its
inherent likelihood-based framework, streamlines statistical inference for
parameter estimation and constraint prediction from 21 cm IM experiments. We
assess the generated HI maps against external simulation data, utilizing
metrics like HI mass probability density functions and power spectra,
validating the prowess of the model across a broad spectrum of mass scales and
spatial dimensions. 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 in Sec. 2.2 and its link to 21 cm physics in Sec. 2.3. Details on
conditional normalizing flows follow in Sec. 2.4. We evaluate the trained
model in Sec. 3.1. 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 Ly${\alpha}$ 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.

##  2 Theoretical Models and Post-Processing

###  2.1 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 $\omega\propto
m^{-1}\lambda_{\text{dB}}^{-2}$, 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. (2022). 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. (2022), we call this proxy classical
FDM (cFDM). For scales corresponding to wavenumbers $k\sim 0.16-80\
h/\text{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\sim 10^{-22}$ eV and halo mass scales around $M\sim 4\times
10^{9}\ M_{\odot}/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_{\text{H}}\sim 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_{\text{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\times 10^{-22},\ 2\times 10^{-21}$ eV. DM halos are identified
using the friends-of-friends (FoF) algorithm with a standard linking length of
$b=0.2\times\text{(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 $\Omega_{m}=0.3089$,
$\Omega_{\Lambda}=0.6911$, $h=H_{0}/100=0.6774$ and $\sigma_{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$.

###  2.2 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,

$u_{\text{eff}}=(1-x)u_{h}+xu_{c},$ (1)

where $u_{\text{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_{\text{eff,TNG}}$ given by

$u_{\text{eff,TNG}}=0.3u_{\text{eff}}+0.7u_{4},$

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

$x=1+\frac{1}{2y}-\sqrt{\frac{1}{y}+\frac{1}{4y^{2}}},$ (2)

where we define

$y=\frac{t_{*}\Lambda_{\text{net}}(\rho,u_{h})}{\rho(\beta
u_{\text{SN}}-(1-\beta)u_{c})}.$ (3)

Here $t_{*}$ is the star formation time scale, $\Lambda_{\text{net}}$ is the
net cooling rate, $\rho$ is the gas density, $u_{\text{SN}}$ is the specific
energy corresponding to a supernova temperature of $T_{\text{SN}}=5.73\times
10^{7}\text{ K}$, and $\beta$ is the mass fraction in massive stars. Note that
we use the IllustrisTNG values of $t_{*}=3.28\text{ Gyr}$ and $\beta=0.226$.
Furthermore, we can convert between specific energies and temperatures using
the equation

$T=\frac{(\gamma-1)\mu}{k_{B}}u,$ (4)

where $\gamma$ is the adiabatic index, $k_{B}$ is the Boltzmann constant, and
$\mu$ is the mean molecular weight given by

$\mu=\frac{4m_{\mathrm{p}}}{1+3X_{H}+4X_{H}(1-x)}.$ (5)

Here, $m_{\mathrm{p}}$ is the mass of a proton and $X_{H}$ is the hydrogen
mass fraction.

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 $\Delta L\approx 2.5\
h^{-1}$Mpc corresponds to a channel width of $\Delta\nu=272$ kHz. Redshift-
space distortions (RSDs) are included by displacing Voronoi cells following
Eq. (12).

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_{\text{H}_{2}}$ in star-forming cells as

$f_{\text{H}_{2}}=\begin{cases}1-\cfrac{0.75s}{1+0.25s}&\text{if $s<2$,}\\\
0&\text{if $s\geq 2$,}\end{cases}$ (6)

where $s$ is given by

$s=\cfrac{\log(1+0.6\chi+0.01\chi^{2})}{0.6\tau_{c}},$ (7)

and

$\chi=0.756(1+3.1Z^{0.365}),\ \ \tau_{c}=\Sigma\sigma_{d}/\mu_{\text{H}}.$ (8)

In the above equations, $Z$ is the gas metallicity (in units of solar
metallicity), $\sigma_{d}=Z\times 10^{-21}\text{ cm}^{2}$ is the cross-section
of dust, $\mu_{\text{H}}=2.3\times 10^{-24}\text{ g}$ is the mean mass per
hydrogen nucleus, and $\Sigma=\rho R$ is the surface density of the gas, where
$\rho$ is the gas density and the volume $V$ of the cell can be translated to
an effective radius $R=(3V/4\pi)^{1/3}$.

Fig. 1 shows 2D HI mass projections of post-processed $\Lambda$CDM and cFDM
simulation data. We choose a frequency channel
$[\nu_{0}-\Delta/2,\nu_{0}+\Delta/2]$ of width $\Delta\nu=272$ kHz around the
central frequency $\nu_{0}=\nu_{21}/(1+z)$, where $\nu_{21}=1420.406$ MHz is
the rest frequency of the 21 cm line. The channel width corresponds to a
comoving width of

$\Delta L=r_{\nu_{21}-\Delta\nu/2}-r_{\nu_{21}+\Delta\nu/2}\approx 2.5\
h^{-1}\text{Mpc},$ (9)

where $r_{\nu}$ is the comoving distance to redshift

$z=\nu_{21}/\nu-1,$ (10)

for the observational frequency $\nu$. We take the same slice of width $\Delta
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).

###  2.3 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 $\bar{x}_{\mathrm{HI}}\approx 10^{-4}$
according to Ly${\alpha}$ and Ly${\beta}$ 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 self-shield against ionizing background
radiation. While most of the HI is found within halos, typically of mass
$M_{\mathrm{h}}\sim 10^{10}-10^{13}\ M_{\odot}$, 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

$\frac{n_{1}}{n_{0}}=3\exp\left(-\frac{h\nu_{21}}{k_{B}T_{s}}\right),$

where $n_{0}$ and $n_{1}$ are the number of ground and excited spin states
respectively and $h$ is Planck’s constant.

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 $\delta T_{\mathrm{b}}\approx 350$ mK are found in DLAs which can
self-shield against external UV radiation.

The differential brightness temperature $T_{b}$ describes the difference
between the CMB temperature $T_{\gamma}$ and the spin temperature $T_{s}$ and
is given by

$\delta
T_{b}(\hat{\mathbf{r}},\nu)=\left[1-e^{-\tau_{21}(\hat{\mathbf{r}},z)}\right]\frac{T_{s}(\hat{\mathbf{r}},z)-T_{\gamma}(z)}{1+z},$

where $\hat{\mathbf{r}}$ is a unit vector that is centred on the observer and
points in the direction of observation on the sky. The observational frequency
$\nu$ is related to the redshift $z$ by the standard relation, Eq. (10). We
have also defined the optical depth $\tau_{21}$ across the 21 cm line at
redshift $z$ (Pritchard & Loeb, 2012; Liu & Shaw, 2020),

$\tau_{21}(\hat{\mathbf{r}},z)=\frac{3hc^{3}A_{10}}{32\pi
k_{B}\nu_{21}^{2}}\frac{x_{\mathrm{HI}}n_{\mathrm{H}}}{(1+z)(dv_{||}/dr_{||})T_{s}},$

where $n_{\mathrm{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\lesssim 100$ K)
or warm ($T\gtrsim 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, low-temperature 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 $\approx 3\times 10^{3}$ K in
the post-reionization era. Due to Ly${\alpha}$ coupling, the spin temperature
of HI is thus much higher than the CMB temperature, and we expect to see the
21 cm line in emission (see e.g. Ando et al., 2021). Therefore, in this
regime, we can assume $T_{\gamma}/T_{s}\approx 0$. Even though certain DLAs at
low redshift $z<1$ can have peak optical depths (as inferred from the 21 cm
line center) of $\tau_{21}\approx 0.5$ or beyond (Curran, 2017; Sadler et al.,
2020), it is customary to Taylor-expand the brightness temperature,

$\delta T_{b}(\hat{\mathbf{r}},\nu)=\frac{3hc^{3}A_{10}}{32\pi
k_{B}\nu_{21}^{2}}\left[\frac{x_{\mathrm{HI}}n_{\mathrm{H}}}{(1+z)^{2}(\mathrm{d}v_{||}/\mathrm{d}r_{||})}\right].$
(11)

In fact, the assumption $T_{\gamma}/T_{s}\approx 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 $\delta 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 $\mathrm{d}v_{||}/\mathrm{d}r_{||}$ 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 $(\mathbf{r})$ to redshift space
$(\mathbf{s})$ through

$\mathbf{s}=\mathbf{r}+\frac{1+z}{H(z)}\mathbf{v}_{||}(\mathbf{r}),$ (12)

where
$H(z)=H_{0}\sqrt{\Omega_{\text{r,0}}(1+z)^{4}+\Omega_{\text{m,0}}(1+z)^{3}+\Omega_{\Lambda,0}}$
is the familiar Hubble parameter and $\mathbf{v}_{||}(\mathbf{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\approx 10$ as a homogeneous, isotropic radiation field,
instead of implementing a more realistic patchy reionization such as the
Sherwood simulations (Puchwein et al., 2023). However, in the post-
reionization era, the assumption of a uniform UV background leads to good
agreement with observations of the mean Ly${\alpha}$ flux, the HI abundance
and the distribution of 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 $\delta T_{\text{b}}$ is
strictly positive throughout the maps, its magnitude spans a large range of
over $7$ orders of magnitude. The maximum brightness temperature reach $\delta
T_{b,\text{max}}\approx 350$ mK and are typically associated with DLAs (see
Sec. 3.4). Note that $\delta T_{b,\text{max}}$ is highly dependent on
numerical resolution. For instance, peak values at $z=4.94$ for the $512^{3}$
runs are notably reduced to $\delta T_{b,\text{max}}\approx 95$ mK.

###  2.4 Conditional Normalizing Flows

The concept of normalizing flows has gained prominence in the field of
probability density estimation, an application first introduced by Dinh et al.
(2015). We now show that the distribution of HI can be effectively learned
using normalizing flows, coupled with rapid sample generation.

Our training dataset consists of pairs $\\{(x_{i},c_{i})\\}_{i=1}^{N}$, where
each $x_{i}$ represents a set of normalized HI maps with dimensions $64\times
64$, and the $c_{i}$ correspond to the logged axion mass parameters
$\log_{10}(m_{a})$ that serve as conditional parameters for the model. By
virtue of the Schrödinger-Vlasov correspondence (Mocz et al., 2018), the
dynamics of FDM and N.B. cFDM converge to CDM dynamics when the de Broglie
wavelength $\lambda_{\text{dB}}$ is significantly smaller than the spatial
resolution $\Delta x$ of our HI maps. We impose

$\frac{\lambda_{\text{dB}}}{\Delta
x}=\frac{h}{m_{a,\text{CDM}}\varv}\frac{1}{\Delta x}<10^{-3}$ (13)

and find $m_{a,\text{CDM}}=10^{-19}$ eV, i.e $-19$ as the CDM conditional
parameter, assuming a local characteristic velocity of DM (in the Madelung
formulation) of $\varv\sim 20$ km/s (Dome et al., 2022). Our full training set
thus contains $40,000$ samples. Our objective now is to approximate the true
underlying probability density function (PDF) $p(x|c)$ using a function
$p_{\boldsymbol{\theta}}(x|c)$.

####  2.4.1 Loss Function

A normalizing flow

$f=f_{1}\circ f_{2}\circ...\circ f_{k}$ (14)

is a composition of invertible functions $(f_{i})_{i=1}^{k}$ that maps a
simple known distribution (such as a spherical multivariate Gaussian noise
$y\sim N(0,\textbf{I})$) to a more complicated distribution (such as the
probability distribution of our HI maps). Note that the functions all
implicitly depend on the parameters $\boldsymbol{\theta}$ of the model. We
define intermediate variables

$h_{i}=f_{1}\circ f_{2}\circ...\circ f_{i}(x,c)$ (15)

with $h_{0}=(x,c)$ and $h_{k}=y$. The key to turning this into an efficient
computational model is to choose the functions $f_{i}$ to be simple with known
analytical inverses and Jacobian determinants. This allows us to write the
loss function of our model as the Kullback-Leibler (KL) divergence between the
underlying PDF $p$ and the model’s estimated PDF $p_{\boldsymbol{\theta}}$:

$L(\boldsymbol{\theta})=D_{KL}(p||p_{\boldsymbol{\theta}}),$ (16)

which we minimize with respect to $\boldsymbol{\theta}$. The KL divergence
between two distributions $p$ and $q$,

$D_{KL}(p||q)=\int_{-\infty}^{\infty}p(x)\log\left(\frac{p(x)}{q(x)}\right)dx,$
(17)

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):

$L(\boldsymbol{\theta})=-\frac{1}{N}\sum_{i=1}^{N}\log
p_{\boldsymbol{\theta}}(x_{i}|c_{i}).$ (18)

Hence our loss function is just the negative log-likelihood. Applying the
standard change of variables formula for PDFs, we get

$\displaystyle\log p_{\boldsymbol{\theta}}(x|c)$ $\displaystyle=\log
p_{\boldsymbol{\theta}}(y)+\log\left|\det\left(\frac{dy}{d(x,c)}\right)\right|$
$\displaystyle=\log
p_{\boldsymbol{\theta}}(y)+\sum_{i=1}^{k}\log\left|\det\left(\frac{dh_{i}}{dh_{i-1}}\right)\right|.$
(19)

In Eq. (19), 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).

Layer | Forward Operation
---|---
Conditional affine injector | $x\mapsto g_{s}(c)\circ x+g_{t}(c),x=(x_{a},x_{b})$
Invertible convolution | $x\mapsto\mathbf{W}x,x=(x_{a},x_{b})$
Conditional affine coupling | $x_{a}\mapsto x_{a},x_{b}\mapsto f_{s}(x_{a},c)\circ x_{b}+f_{t}(x_{a},c)$
Table 1: 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., 2021). Note how the model is made
conditional by introducing the conditional parameter $c$ into the affine
injector and affine coupling layers.

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).$ (20)

####  2.4.2 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., 2017) and NICE (Dinh et al.,
2015). 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\times 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., 2017) 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 to incorporate conditioning on
parameters. Table 1 shows the three key layers that make up a conditional flow
block.

Hyperparameter | Value
---|---
Flow blocks | 6
Flow steps per block | 12
Hidden width | 20
Learning rate | $3\times 10^{-4}$
Sigmoid rescaling $s$ | $0.5$
Epochs | 1000
Table 2: Hyperparameter values used for our model. Both the conditional affine
coupling as well as the affine injector have one hidden layer of width given
above. The sigmoid rescaling parameter $s$ refers to the constant in Eq. (23).

The inverse flow and thus sample generation can be prone to numerical
instabilities (Behrmann et al., 2021), even if the loss curve plateaus as
expected and training is successful. To stabilize the inverses against small
numerical errors, we choose small Lipschitz constants $L_{i}$ for the
functions that constitute the flow and their respective inverses $g_{i}$. By
rearranging the definition of Lipschitz continuity,

$\frac{\norm{g_{i}(x_{1})-g_{i}(x_{2})}}{\norm{x_{1}-x_{2}}}\leq L_{i},$ (21)

a smaller Lipschitz constant $L_{i}$ thus implies that $g_{i}$ has small
gradients (if differentiable). An upper bound on the Lipschitz constant for
the inverse of an affine coupling block (see Table 1) is established in
Behrmann et al. (2021). The bound depends linearly on the reciprocal of the
sigmoid function and its derivative. However, the reciprocal of the sigmoid
function grows exponentially to $\pm$infinity for large inputs close to $0,1$,
respectively. Thus, the Lipschitz constant becomes potentially unbounded, and
so do the numerical errors associated with this layer. The sigmoid function

$S(x)=\frac{1}{1+\exp(-x)},$ (22)

appears as the function $f_{s}$ in the affine coupling layer of Glow (see
Table 1). We thus modify the affine coupling layer to use a linearly rescaled
sigmoid function,

$\tilde{S}(x)=(1-s)S(x)+s,$ (23)

where $s$ is a constant we set. $\tilde{S}$ takes values on the interval
$(s,1)$ rather than $(0,1)$ and therefore assures that the relevant Lipschitz
constant is bounded from above for large enough values of $s>0$.

The key hyperparameters used for training are shown in Table 2. To allow the
maximum amount of splitting (Dinh et al., 2017) of the $64^{2}$ images
throughout the flow, we employ $6$ flow blocks. This enables effective
utilization of spatial structures at various scales during training. $12$
flows per block and a hidden layer width of $20$ make sure the model has
enough parameters to learn the key features without over- or underfitting. On
average, training took $6$ hours per $100$ epochs on an NVIDIA A100 GPU. The
final model that is presented in this work was trained for $1000$ epochs.

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\times 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.

####  2.4.3 Data Pre-Processing

We pick $10,000$ random $64^{3}$ subcubes from the full 3D HI distribution in
each of the four cosmologies. The subcubes are then projected to get $64^{2}$
HI maps. These $64^{2}$ maps have a comoving length of $2.5\ h^{-1}$Mpc; this
length scale is chosen to balance the requirement of a large number of samples
while keeping the individual maps sufficiently large to capture DM imprints.
Using projected maps leads to more stable and less computationally expensive
training.

Figure 4: Statistical properties of simulation data (dotted lines) and HIGlow-
generated 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\times 10^{-22}$
eV (see Fig. 3). This case does not have corresponding simulation data. The
Nyquist frequency $k_{\text{Ny}}=\pi N^{1/3}/L_{\text{box}}$ is shown as a
vertical line.

If $X$ denotes the unnormalized map from our training data, we first take the
logarithm of each map, $x=\log_{10}(X)$. Now we sigmoid-normalize our data
onto the interval $(-1,1)$ by

$\tilde{x}=\frac{2}{1+\exp(-\hat{x})}-1,$ (24)

where $\bar{x}$ and $\sigma_{x}$ is the mean and standard deviation of the
whole training dataset, respectively, and $\hat{x}=(x-\bar{x})/\sigma_{x}$.
Fig. 3 shows some sigmoid-normalized $64\times 64$ HI training maps in the top
four rows across the CDM and cFDM cosmologies at $z=4.94$. The bottom four
rows show generated HIGlow HI samples with identical seeds across the
cosmologies. There are no diverging pixels or other visual anomalies. Without
the rescaled sigmoid in the conditional affine coupling layer, Eq. (23), an
alternative approach would have been to apply image clipping, constraining
pixel values to fall within predefined maximum and minimum ranges, after each
of the six flow blocks in the generative reverse flow. However, this would
have come at the cost of losing important image features with each subsequent
clipping operation.

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\times 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\times 10^{-22}$ eV
samples, aligning with our expectations. We have conducted comprehensive
validation tests, the results of which are outlined in the following section.

##  3 Post-Reionization HI Distributions

###  3.1 HIGlow Modeling of HI

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_{\text{HI}}\approx
10^{5}\ M_{\odot}/h$. In the range $m_{\text{HI}}=10^{5}-10^{7}\ M_{\odot}/h$,
HIGlow captures the distribution less accurately, which is a result of the
rarity of such high-mass pixels. For instance, $m_{\text{HI}}=10^{6}\
M_{\odot}/h$ pixels are $6-7$ orders of magnitude less likely than pixels with
$m_{\text{HI,peak}}\approx 20\ M_{\odot}/h$. The PDF bump around this peak
value broadly corresponds to the Ly${\alpha}$ forest in the post-reionization
era (also see Zamudio-Fernandez et al., 2019).

Recall that we sigmoid-normalise the data before training HIGlow following Eq.
(24). The power spectra of Fig. 4 are also shown in sigmoid space since the
trends with axion mass $m$ are more clear. 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-sigmoid-normalized 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\times 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\times 10^{-22}$ eV, as 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\times 10^{-22}$
eV, the generated distribution effectively showcases the interpolation prowess
in latent space. While a direct comparison to actual simulations at $m=3\times
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.

###  3.2 Overall HI Abundance

Here we investigate the overall HI abundance
$\Omega_{\text{HI}}=\bar{\rho}_{\text{HI}}(z)/\rho_{\mathrm{crit}}^{0}$
relative to the critical density of the Universe today
$\rho_{\mathrm{crit}}^{0}$, with $\bar{\rho}_{\text{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\approx 0$) HI at around
$\Omega_{\mathrm{HI}}\sim(3.9\pm 0.6)\times 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\sim 0.8$ to be
$\Omega_{\mathrm{HI}}b_{\mathrm{HI}}\sim 6.2^{+2.3}_{-1.5}\times 10^{-4}$
(Switzer et al., 2013). More indirect measurements of DLAs allow us to
constrain the HI abundance up to $z\sim 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}\propto\Omega_{\text{HI}}^{2}$.

We show the value of $\Omega_{\text{HI}}(z)$ in the four simulated cosmologies
in Fig. 5 in the redshift range $z=3.42-4.94$. 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

$f_{\text{HI}}(N_{\text{HI}},X)=\frac{\mathrm{d}^{2}N(N_{\text{HI}})}{\mathrm{d}N_{\text{HI}}\mathrm{d}X},$
(25)

where $N$ is the number of lines of sight with column densities between
$N_{\text{HI}}$ and $N_{\text{HI}}+\mathrm{N_{\text{HI}}}$, and

$\mathrm{d}X=\frac{H_{0}}{H(z)}(1+z)^{2}\mathrm{d}z$ (26)

is the so-called absorption distance. The absorption distance has the property
that absorbers with non-evolving number density $n_{a}$ and cross-section
$\sigma_{a}$ (both in proper length units) will produce a constant number of
absorption lines per unit absorption distance, i.e.
$\mathrm{d}N/\mathrm{d}X=\text{const}$.

Figure 5: Comparison of mock and observed cosmic HI density estimates. The
abundance parameter
$\Omega_{\text{HI}}=\bar{\rho}_{\text{HI}}(z)/\rho_{c}^{0}$ 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\sigma$ 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: 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 per Marsh
(2016) and the Nyquist frequency $k_{\text{Ny}}=\pi N^{1/3}/L_{\text{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_{\text{HI}}^{2}(k)=P_{\text{HI}}(k)/P_{\text{tot}}(k)$ increase
for decreasing axion mass $m$.

The distribution $f_{\text{HI}}(N_{\text{HI}},X)$ spans ten orders of
magnitude from around $N_{\text{HI}}=10^{12}$ cm-2, below which Ly${\alpha}$
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 $\simeq 10^{17.3}$ cm-2 are called the Ly${\alpha}$ forest,
and physically represent the low-density and highly ionized gas that resides
in the IGM. Systems with column densities above $N_{\text{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${\alpha}$ forest
that accounts for most of the gas (neutral + ionized).

The total gas density can be inferred from the HI CDDF via

$\Omega_{\mathrm{g}}^{\mathrm{DLA}}(X)\mathrm{d}X=\frac{\mu
m_{\mathrm{p}}H_{0}}{c\rho_{\mathrm{crit}}^{0}}\int_{10^{20.3}}^{\infty}N_{\mathrm{HI}}f_{\mathrm{HI}}(N_{\mathrm{HI}},X)\mathrm{d}X,$
(27)

where $\mu=1.3$ is the mean molecular mass of the gas. In the discrete limit,
$\Omega_{\mathrm{g}}^{\mathrm{DLA}}$ is given by

$\Omega_{\mathrm{g}}^{\mathrm{DLA}}=\frac{\mu
m_{\mathrm{p}}H_{0}}{c\rho_{\mathrm{crit}}^{0}}\frac{\sum
N_{\mathrm{HI}}}{\Delta X},$ (28)

where the sum is calculated for systems with $\log_{10}N_{\mathrm{HI}}\geq
20.3$ along lines of sight with a total pathlength of $\Delta X$. As noted in
(Zafar et al., 2013; Crighton et al., 2015), we can convert the gas mass in
DLAs $\Omega_{\mathrm{g}}^{\mathrm{DLA}}$ to the neutral hydrogen mass via

$\Omega_{\mathrm{HI}}=\delta_{\mathrm{HI}}\Omega_{\mathrm{g}}^{\mathrm{DLA}}/\mu,$
(29)

where $\mu$ accounts for the mass of helium and $\delta_{\mathrm{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 $\Delta X$ in Eq. (27) 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)111The thermal relic WDM masses explored
in Carucci et al. (2015) correspond to axion masses of $m\in[1.6\times
10^{-22},5.5\times 10^{-21}]$ eV..

We note that similar to Villaescusa-Navarro (2018), the overall HI abundance
depends on resolution. The lower-resolution $N=512^{3}$ simulations contain
$\sim 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}\sim 10^{9}\ M_{\odot}/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.

###  3.3 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,

$P_{X}=\frac{1}{N_{\text{modes}}}\sum_{\mathbf{k}\in
k}\delta_{X}(\mathbf{k})\delta_{X}^{\ast}(\mathbf{k}),$ (30)

with $X=\\{\text{tot},\text{HI}\\}$ and $N_{\text{modes}}$ being the number of
modes lying in the spherical shell $k$ of width $\delta
k=2\pi/L_{\text{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_{\text{Ny}}=\pi N^{1/3}/L_{\text{box}}$, where
$N^{1/3}$ is the number of cells per box side.

Figure 7: 21 cm power spectrum in CDM and cFDM cosmologies at $z=4.94$. We
show results for the ‘dimensionless’ power spectrum
$\Delta^{2}_{21}(k)=k^{3}P_{21}(k)/(2\pi^{2})$, in both real space (dotted
lines) and redshift space (solid lines). Black crosses denote the expected
SKA1-Low $1\sigma$ system temperature noise $\Delta^{2}_{\text{noise}}$ for
$1080$ hours of observation in a $8$ MHz bandwidth with $290$ frequency
channels as per Table 5, for an optimistic foreground wedge model. The
vertical dashed line denotes the Nyquist frequency $k_{\text{Ny}}=\pi
N^{1/3}/L_{\text{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
$\Delta^{2}_{\text{noise}}/\Delta^{2}_{21,\text{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$.

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 small-scale 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_{\text{HI}}(k)=\sqrt{\frac{P_{\text{HI}}(k)}{P_{\text{tot}}(k)}}$ (31)

describes the bias between the distributions of HI and the total matter. While
the HI bias defined using Eq. (31) suffers from stochasticity (shot noise), it
is closer to observations than results inferred using the alternative
definition involving the HI-matter cross-correlation,
$b_{\text{HI}}(k)=P_{\text{HI-tot}}(k)/P_{\text{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_{\text{HI}}\approx 2.5$
in redshift space and $b_{\text{HI}}\approx 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_{\text{HI}}\approx 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 $\approx
2.5$.

The clustering of the 21 cm signal which is modeled following Sec. 2.3 is
shown in Fig. 7. Its magnitude $\propto\Omega_{\text{HI}}^{2}(z)$ in the post-
reionization 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 is shown for instrument parameters from Table 5,
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\approx 9$ $h$/Mpc in CDM and the
less extreme cFDM models $m\sim 10^{-21}$ eV, while it can be detected up to
$k\approx 15$ $h$/Mpc in the more extreme $m=10^{-22}$ eV model.

Since we assume that $T_{\gamma}/T_{s}\approx 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 $\Delta^{2}_{\text{noise}}/\Delta^{2}_{21,\text{CDM}}$
for SKA1-Low, representing the $1\sigma$ 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\sigma$
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\sim 10^{-21}$ eV with $1080$ hours of observations on
small scales $k\approx 70$ $h$/Mpc.

###  3.4 HI Column Density Distributions

Another quantity commonly employed to study the distribution of HI in the
post-reionization era is the HI CDDF, Eq. (25). As mentioned in Sec. 2.2, HI
column densities of Ly${\alpha}$ clouds, LSS, sub-DLAs and DLAs are inferred
observationally from quasar spectra.

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\sigma$ 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.

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_{\text{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 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_{\text{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:

$N_{\text{HI},i}=2\frac{m_{\text{HI}}}{m_{\text{H}}}\int_{0}^{l_{\text{max}}}W(r,h_{i})\mathrm{d}l,$
(32)

where $N_{\text{HI},i}$ is the column density due to particle $i$, having HI
mass $m_{\text{HI}}$ and SPH smoothing length $h_{i}$ while $m_{\text{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_{\text{max}}^{2}=h_{i}^{2}-b^{2}$. Our regular grid consists of
$10,000\times 10,000$ points, i.e. the number of LOSs is $10^{8}$. In view of
our small box size $L_{\text{box}}=40\ h^{-1}$Mpc, the probability of
encountering more than a single absorber with a large column density $\sim
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 $\sim 10^{19}$ cm-2.

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\times
10^{-21}$ eV (top right) and $m=7\times 10^{-22}$ eV (bottom left)
cosmologies, we fit a power law + exponential function (34) while for the
$m=10^{-22}$ eV (bottom right) cosmology we fit a double power law, Eq. (35).
The vertical dashed line denotes the half-mode mass $M_{1/2}$ (Marsh, 2016).
The small secondary peak in the $m=7\times 10^{-22}$ eV model morphs into a
bimodal distribution in the $m=10^{-22}$ eV cosmology.

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) 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 (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.

For CDM and the less extreme cFDM models $m=7\times 10^{-22}$ eV and
$m=2\times 10^{-21}$ eV, we find good agreement with the observations. The
strongest DLA absorbers are not modeled well in any of the simulated
cosmologies, which can be traced to our moderate box size of
$L_{\text{box}}=40\ h^{-1}$Mpc. The more extreme $m=10^{-22}$ eV model is
visually at odds with the measurements since the predicted CDDF undercuts
observational estimates for column densities above $N_{\text{HI}}=2\times
10^{21}$ cm-2. With larger data sets and better mean flux measurements, one
might be able to put competitive astrophysical constraints on FDM models.

Table 3: Best-fit parameters $M_{0}$ in Eq. (34) for CDM, $m=2\times 10^{-21}$ eV and $m=7\times 10^{-22}$ eV cosmologies and $M_{\mathrm{break}}$, $\beta$ in Eq. (35) for the $m=10^{-22}$ eV cosmology across redshifts $z=3.42-4.94$. |  $M_{0}\ [M_{\odot}/h]$ |  $M_{\mathrm{break}}\ [M_{\odot}/h]$ |  $\beta$
---|---|---|---
$z=3.42$ |  $3.6\times 10^{9}$ |  $1.5\times 10^{9}$ |  $0.52$
$z=3.88$ |  $2.6\times 10^{9}$ |  $1.6\times 10^{9}$ |  $0.51$
$z=4.38$ |  $1.9\times 10^{9}$ |  $1.7\times 10^{9}$ |  $0.48$
$z=4.94$ |  $1.6\times 10^{9}$ |  $1.6\times 10^{9}$ |  $0.45$

###  3.5 DLA Cross-Sections

Next, we focus on absorbers with high column densities above
$N_{\text{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

$\sigma_{\text{DLA}}=\pi
R_{\text{vir}}^{2}\left(\frac{n_{\text{DLA}}}{n_{\text{tot}}}\right),$ (33)

where $R_{\text{vir}}$ is the halo virial radius, $n_{\text{DLA}}$ is the
number of LOSs with column density above $10^{20.3}$ cm-2 and $n_{\text{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_{\text{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.

Table 4: DLA bias $b_{\text{DLA}}$, estimated using Eq. (36), in CDM and cFDM cosmologies across redshifts $z=3.42-4.94$. |  $\mathrm{CDM}$ |  $2\times 10^{-21}\ \mathrm{eV}$ |  $7\times 10^{-22}\ \mathrm{eV}$ |  $10^{-22}\ \mathrm{eV}$
---|---|---|---|---
$z=3.42$ |  $1.8$ |  $2.0$ |  $2.1$ |  $2.6$
$z=3.88$ |  $1.9$ |  $2.2$ |  $2.3$ |  $2.9$
$z=4.38$ |  $2.1$ |  $2.3$ |  $2.6$ |  $3.3$
$z=4.94$ |  $2.3$ |  $2.6$ |  $2.9$ |  $3.8$

We show results in Fig. 9. We find that the mean DLA cross-section increases
with halo mass, which apart from the $m=10^{-22}$ eV cosmology is well fitted
by the following function

$\sigma(M|N_{\text{HI,DLA}},z)=A\left(\frac{M}{h^{-1}M_{\odot}}\right)^{\alpha}\left(1-e^{-(M/M_{0})^{\beta}}\right).$
(34)

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 $\alpha=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, $A\cdot M_{0}=14100\
h^{-3}\text{kpc}^{2}M_{\odot}$. The cutoff exponent is well approximated by
$\beta=0.85\cdot\log_{10}(N_{\text{HI,DLA}}/\text{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\times 10^{9}\
h^{-1}M_{\odot}$ in all three cosmologies. In Table 3, 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\times 10^{-22}$ eV model, the $m=10^{-22}$ eV cosmology exhibits a
bimodal distribution at $z=4.94$. The median cross-section is significantly
higher than in the other three cosmologies for halo masses below $M\approx
3\times 10^{9}\ M_{\odot}/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\approx 5$ is around $30$% higher in the
$m=10^{-22}$ eV cosmology than for CDM (Dome et al., 2023). The power law +
cutoff parametrization of Eq. (34) thus gives a poor fit for the $m=10^{-22}$
eV model, and instead we adopt a double power law

$\sigma(M|N_{\text{HI,DLA}},z)=\begin{cases}&A\left(\frac{M}{h^{-1}M_{\odot}}\right)^{\alpha}\
\ \text{if}\ \ M>M_{\mathrm{break}}\\\
&B\left(\frac{M}{h^{-1}M_{\odot}}\right)^{\beta}\ \ \text{else}.\end{cases}$
(35)

The dependence of $\sigma_{\mathrm{DLA}}$ on cosmology at the high-mass end
(i.e. $M>M_{\mathrm{break}}$) is so weak that the normalization $A$ can still
be well approximated by $A\cdot M_{0}=14100\ h^{-3}\text{kpc}^{2}M_{\odot}$,
where $M_{0}$ is taken from Eq. (34). The normalization at the low-mass end
$B$ is determined by ensuring that the two power laws intersect at
$M_{\mathrm{break}}$. We provide the best-fit values for $M_{\mathrm{break}}$
and $\beta$ in Table 3.

Using the mean DLA cross-section, we can estimate the DLA bias as

$b_{\text{DLA}}(z|N_{\text{HI,DLA}})=\frac{\int_{0}^{\infty}b(M,z)n(M,z)\sigma(M|N_{\text{HI,DLA}},z)\mathrm{d}M}{\int_{0}^{\infty}n(M,z)\sigma(M|N_{\text{HI,DLA}},z)\mathrm{d}M},$
(36)

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 $\nu=\delta_{c}/\sigma(M)$, where
$\delta_{c}=1.686$ is the spherical collapse threshold and $\sigma^{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 $\sigma^{2}(M)$ in
the respective cosmology (Dunstan et al., 2011).

In Table 4, 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_{\text{HI}}>N_{\text{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_{\text{DLA}}\simeq
b_{\text{HI}}$, following the theoretical arguments in Castorina &
Villaescusa-Navarro (2017) and comparing the values in Table 4 to the
$b_{\text{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_{\text{DLA}}=2.3$ for CDM to
$b_{\text{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${\alpha}$ forest, Pérez-Ràfols et al.
(2018) measured a DLA bias factor of $b_{\text{DLA}}=1.99\pm 0.11$, while
using cross-correlations with CMB lensing data Alonso et al. (2018) and Lin et
al. (2020) measured $b_{\text{DLA}}=2.6\pm 0.9$ and
$b_{\text{DLA}}=1.37^{+1.30}_{-0.92}$, respectively. All three estimates are
for a median DLA redshift of $z_{\mathrm{median}}=2.3$ and agree well with our
CDM estimates from Table 4, 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. (36), the observed DLA bias could be used to put
constraints on FDM-like cosmologies.

##  4 Mock Radio Maps for SKA-Low

Configuration |
---|---
Array design | 512 compact core
Station diameter $D$ [m] | 35
Integration time $t_{\mathrm{int}}$ [s] | 60
Bandwidth $B$ [MHz] | 8
Number of channels $N_{\text{chan}}$ | 29/290a
Frequency resolution $\Delta\nu$ [kHz] | 272
Number of observing hours/days | 6/180
Redshift/frequency [MHz] | 4.94/238.9
Angular resolution $\theta_{\mathrm{A}}$ | 0.5’, 1’, 2’
Primary beam’s FWHM $\theta_{\mathrm{max}}$ | 3.12($\nu/200$ MHz)-1 deg
Single baseline noise error $\sqrt{\langle|N|^{2}\rangle}$ [Jy]/[mK] | 0.55/134
aFor imaging: $N_{\text{chan}}=29$; for 21 cm power spectra:
$N_{\text{chan}}=290$.
Table 5: 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. (2020). Natural sensitivities $A_{\mathrm{eff}}/T_{\mathrm{sys}}$ as a function of frequency are tabulated in Braun et al. (2019). Table 6: Prospects of detecting the brightest HI peaks using SKA1-Low at $z=4.94$. We provide the peak flux $I_{\mathrm{max}}$, noise rms $\sqrt{\langle|N|^{2}\rangle}$ (both in $\mu$Jy/beam) and SNR values in CDM and cFDM cosmologies for channel width $\Delta\nu=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\times 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. | CDM |  $m=2\times 10^{-21}$ eV |  $m=7\times 10^{-22}$ eV |  $m=3\times 10^{-22}$ eV |  $m=10^{-22}$ eV |
---|---|---|---|---|---|---
$\theta_{A}$ |  $I_{\mathrm{max}}$ |  $\mathrm{SNR}$ |  $I_{\mathrm{max}}$ |  $\mathrm{SNR}$ |  $I_{\mathrm{max}}$ |  $\mathrm{SNR}$ |  $I_{\mathrm{max}}$ |  $\mathrm{SNR}$ |  $I_{\mathrm{max}}$ |  $\mathrm{SNR}$ |  $\sqrt{\langle|N|^{2}\rangle}$
$0.5^{\prime}$ |  $0.82$ |  $0.88$ |  $0.78$ |  $0.84$ |  $0.76$ |  $0.82$ |  $0.74$ |  $0.80$ |  $0.70$ |  $0.75$ |  $0.93$
$1.0^{\prime}$ |  $1.39$ |  $1.81$ |  $1.31$ |  $1.70$ |  $1.27$ |  $1.65$ |  $1.23$ |  $1.60$ |  $1.17$ |  $1.52$ |  $0.77$
$2.0^{\prime}$ |  $2.08$ |  $3.25$ |  $1.97$ |  $3.08$ |  $1.93$ |  $2.84$ |  $1.82$ |  $3.02$ |  $1.70$ |  $2.66$ |  $0.64$

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
$\Delta\nu=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

$I_{\nu}(\hat{\mathbf{r}},\nu)=\frac{2\nu^{2}}{c^{2}}k_{\text{B}}\delta
T_{\text{b}}(\hat{\mathbf{r}},\nu),$ (37)

and project the $I_{\nu}(\hat{\mathbf{r}},\nu)$ field onto a regular grid of
resolution $1024\times 1024$ in the XY plane.

###  4.1 Angular Resolution

The intrinsic resolution of our simulation $\Delta x=L_{\mathrm{box}}/N\approx
40\ h^{-1}$kpc corresponds to angular scales

$\Delta\theta=\frac{\Delta x}{D_{c}(z)}\approx 1.5^{\prime\prime},$ (38)

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_{\mathrm{max}}$ via
$\theta_{A}=\lambda/B_{\mathrm{max}}$, where $\lambda$ is the signal
wavelength at the observation redshift $z$. Since the maximum baseline of
SKA1-Low is $B_{\mathrm{max}}=65.4$ km corresponding to $\approx
4.0^{\prime\prime}$, SKA1-Low could get close to such angular resolutions.
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
$\theta_{A}=0.5^{\prime},1^{\prime},2^{\prime}$ corresponding to baselines
$B_{\mathrm{max}}=8.6,4.3,2.2$ km.

We account for the chosen instrumental resolution by smoothing the $I_{\nu}$
field with a Gaussian window of angular radius $\theta_{A}$, with $\theta_{A}$
being the FWHM of the synthesized beam given in Table 5. 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_{\nu}\rightarrow I_{\nu}-\langle
I_{\nu}\rangle$. Finally, we multiply the $I_{\nu}$ map with the beam solid
angle

$\Delta\Omega=\left(\frac{\pi}{4\ln 2}\right)\theta_{A}^{2}$ (39)

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_{\nu}$, and report results in Table 6.

###  4.2 Thermal Noise

Figure 10: Mock radio maps for CDM with a frequency channel of width
$\Delta\nu=272$ kHz and synthesized beam FWHM $\theta_{A}=1^{\prime}$ (top
row) and $\theta_{A}=2^{\prime}$ (bottom row). Left column: Specific intensity
field $I_{\nu}(\hat{\mathbf{r}})$ for $N_{\text{pix}}=1024$ corresponding to
our simulation resolution. Right column: Specific intensity field overlaid
with instrumental noise for $N_{\text{pix}}=12$ and $N_{\text{pix}}=25$,
respectively. We construct the noise map following a four-step procedure
outlined in Sec. 4.2.

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)

$\sqrt{\langle|N|^{2}\rangle}\mathrm{[Jy]}=\frac{2k_{\mathrm{B}}T_{\mathrm{sys}}}{A_{\mathrm{eff}}\sqrt{\Delta\nu
t_{\mathrm{int}}}},$ (40)

where $t_{\mathrm{int}}$ is the integration time to observe a single
visibility at frequency resolution $\Delta\nu$, $A_{\mathrm{eff}}$ is the
effective area of each station and $T_{\mathrm{sys}}$ is the system
temperature which combines the receiver temperature and the sky temperature.
Their combination $A_{\mathrm{eff}}/T_{\mathrm{sys}}$, the natural
sensitivity, is frequency dependent and we retrieve the most up-to-date tables
from Braun et al. (2019). Eq. (40) provides the flux error for a single
baseline which we convert to Kelvin using

$\sqrt{\langle|N|^{2}\rangle}\mathrm{[K]}=\frac{\lambda^{2}}{2k_{\mathrm{B}}\Omega}\times\sqrt{\langle|N|^{2}\rangle}\mathrm{[Jy]},$
(41)

where $\Omega=\pi/(4\ln 2)\theta_{\mathrm{max}}^{2}$ is the full beam area. We
present rms noise errors for a single baseline in Table 5, along with the
SKA1-Low instrument parameters.

To generate a 2D thermal noise realization we follow four steps:

  1. 1.

We first generate Gaussian random noise (both the real and imaginary parts)
with zero mean and rms $\sqrt{\langle|N|^{2}\rangle}$ from Eq. (40) in the
$N_{\mathrm{pix}}\times N_{\mathrm{pix}}$ grid in Fourier space, where
$N_{\mathrm{pix}}=\theta_{\mathrm{max}}/\theta_{A}$.

  2. 2.

We compute the daily $uv$ coverage $N_{uv}$ which represents the total number
of baselines that observe a given $uv$ pixel using 21cmSense Pober et al.
(2013, 2014). 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/\sqrt{N_{uv}^{ij}}$ at non-zero $uv$ coverage pixels.

  3. 3.

By averaging over the total number of days observed, we reduce the noise
further by a factor of $1/\sqrt{N_{\mathrm{days}}}$.

  4. 4.

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
$\theta_{A}$ obtained using the above four-step procedure are given in Table
6. 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 $\theta_{A}=2^{\prime}$ is modest (SNR
$>3$), for $\theta_{A}=1^{\prime}$ 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 5 are moderate at
$\theta_{A}=2^{\prime}$, with a rapidly declining SNR for lower values of
$\theta_{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., 2019, 2020), 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,
$\theta_{A}=1^{\prime}$ and $\theta_{A}=2^{\prime}$. Notably, a prominent peak
(SNR$=3.25$) can be observed in the lower half for $\theta_{A}=2^{\prime}$,
corresponding to fluxes per beam around $2\mu$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
$\theta_{A}=1^{\prime}$, the peak becomes more challenging to distinguish
(SNR=$1.81$).

###  4.3 HIGlow Interpolation

The imaging prospects for cFDM cosmologies are similar to CDM, see peak flux
and SNR values in Table 6. At our fiducial redshift of $z=4.94$, we find a
trend of decreasing peak flux $I_{\mathrm{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 demonstrate 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\times 10^{-22}$ eV (see Sec. 3.1) and quote the maximum peak flux
$I_{\text{max}}$ among all samples in Table 6. We find that the values are in
between the $m=7\times 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.

##  5 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 non-cold DM fraction (see e.g. Giri & Schneider,
2022). In this work, we model the distribution of neutral hydrogen (HI) using
high-resolution hydrodynamical $N$-body simulations in CDM and three instances
of cFDM cosmologies with axion masses $m=10^{-22},7\times 10^{-22},2\times
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:

  1. (a)

We observe a trend of decreasing HI abundance $\Omega_{\text{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\sim 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 $\Omega_{\text{HI}}$
directly.

  2. (b)

The HI bias $b_{\text{HI}}^{2}(k)=P_{\text{HI}}(k)/P_{\text{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\sigma$ confidence level on all scales $k<80\ h/$Mpc resolved by
the simulations.

  3. (c)

We show that the HI frequency distribution $f_{\text{HI}}(N_{\text{HI}},X)$
decreases in cFDM cosmologies compared to CDM. Its weak dependence on redshift
in the post-reionization era allows to compare to a range of DLA and sub-DLA
observations covering a wide redshift range. They show very good agreement
with predictions from CDM and the less extreme cFDM models $m\sim 10^{-21}$
eV, with tensions appearing for the $m=10^{-22}$ eV model.

  4. (d)

DLA cross-sections exhibit a simple power law + low-mass exponential cutoff
trend for cFDM models $m\sim 10^{-21}$ eV that is well known for CDM models
(Villaescusa-Navarro et al., 2018). For more extreme models with $m\sim
10^{-22}$ eV, we show that a double power law parametrization provides a
better fit. The high median cross-section 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_{\text{DLA}}$ and recover the
well-known result of increasing DLA bias with increasing redshift. In
addition, we find the redshift-independent trend that the DLA bias increases
with decreasing axion mass $m$, from $b_{\text{DLA}}=2.3$ for CDM to
$b_{\text{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.

  5. (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
$\theta_{A}\approx 2$, with a rapidly declining SNR for lower values of
$\theta_{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$.

  6. (f)

We demonstrate the ability of normalizing flows, a generative ML model
introduced by Agnelli et al. (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 a 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\times 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 post-reionization 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. Advances
in cosmological simulations 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.

##  6 Acknowledgements

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 addressing pixel divergence issues and 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.

##  7 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.

## References

  * Agnelli et al. (2010) Agnelli J. P., Cadeiras M., Tabak E. G., Turner C. V., Vanden-Eijnden E., 2010, Multiscale Modeling & Simulation, 8, 1784
  * Alonso et al. (2018) Alonso D., Colosimo J., Font-Ribera A., Slosar A., 2018, J. Cosmology Astropart. Phys., 2018, 053
  * Ando et al. (2021) Ando R., Nishizawa A. J., Shimizu I., Nagamine K., 2021, MNRAS, 507, 2937
  * Arvanitaki et al. (2010) Arvanitaki A., Dimopoulos S., Dubovsky S., Kaloper N., March-Russell J., 2010, Phys. Rev. D, 81, 123530
  * Bandura et al. (2019) Bandura K., et al., 2019, Packed Ultra-wideband Mapping Array (PUMA): A Radio Telescope for Cosmology and Transients (arXiv:1907.12559)
  * Barkana (2018) Barkana R., 2018, Nature, 555, 71
  * Barnes et al. (2001) Barnes D. G., et al., 2001, MNRAS, 322, 486
  * Bauer et al. (2021) Bauer J. B., Marsh D. J. E., Hložek R., Padmanabhan H., Laguë A., 2021, MNRAS, 500, 3162
  * Behrmann et al. (2021) Behrmann J., Vicol P., Wang K.-C., Grosse R., Jacobsen J.-H., 2021, Understanding and Mitigating Exploding Inverses in Invertible Neural Networks (arXiv:2006.09347)
  * Bevins et al. (2022) Bevins H., Handley W., Lemos P., Sims P., de Lera Acedo E., Fialkov A., 2022, arXiv e-prints, p. arXiv:2207.11457
  * Bird et al. (2017) Bird S., Garnett R., Ho S., 2017, MNRAS, 466, 2111
  * Bosman et al. (2021) Bosman S. E. I., Ďurovčíková D., Davies F. B., Eilers A.-C., 2021, MNRAS, 503, 2077
  * Boyarsky et al. (2019) Boyarsky A., Iakubovskyi D., Ruchayskiy O., Rudakovskyi A., Valkenburg W., 2019, Phys. Rev. D, 100, 123005
  * Braun et al. (2019) Braun R., Bonaldi A., Bourke T., Keane E., Wagg J., 2019, arXiv e-prints, p. arXiv:1912.12699
  * Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495
  * Bullock & Boylan-Kolchin (2017) Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343
  * Carucci (2018) Carucci I. P., 2018, in Journal of Physics Conference Series. p. 012003 (arXiv:1712.04022), doi:10.1088/1742-6596/956/1/012003
  * Carucci et al. (2015) Carucci I. P., Villaescusa-Navarro F., Viel M., Lapi A., 2015, J. Cosmology Astropart. Phys., 2015, 047
  * Castellano et al. (2023) Castellano M., Menci N., Romanello M., 2023, arXiv e-prints, p. arXiv:2301.03854
  * Castorina & Villaescusa-Navarro (2017) Castorina E., Villaescusa-Navarro F., 2017, MNRAS, 471, 1788
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chang et al. (2008) Chang T.-C., Pen U.-L., Peterson J. B., McDonald P., 2008, Phys. Rev. Lett., 100, 091303
  * Chua et al. (2019) Chua K. T. E., Pillepich A., Vogelsberger M., Hernquist L., 2019, MNRAS, 484, 476
  * Collaboration et al. (2022) Collaboration C., et al., 2022, Detection of Cosmological 21 cm Emission with the Canadian Hydrogen Intensity Mapping Experiment (arXiv:2202.01242)
  * Combes (2021) Combes F., 2021, in Siebert A., et al., eds, SF2A-2021: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. pp 197–202 (arXiv:2107.03915), doi:10.48550/arXiv.2107.03915
  * Corasaniti et al. (2017) Corasaniti P. S., Agarwal S., Marsh D. J. E., Das S., 2017, Phys. Rev. D, 95, 083512
  * Crighton et al. (2015) Crighton N. H. M., et al., 2015, MNRAS, 452, 217
  * Curran (2017) Curran S. J., 2017, A&A, 606, A56
  * Del Popolo & Le Delliou (2017) Del Popolo A., Le Delliou M., 2017, Galaxies, 5, 17
  * Demirtas et al. (2018) Demirtas M., Long C., McAllister L., Stillman M., 2018, arXiv e-prints, p. arXiv:1808.01282
  * Diemer et al. (2018) Diemer B., et al., 2018, ApJS, 238, 33
  * Diemer et al. (2019) Diemer B., et al., 2019, MNRAS, 487, 1529
  * Dinh et al. (2015) Dinh L., Krueger D., Bengio Y., 2015, NICE: Non-linear Independent Components Estimation (arXiv:1410.8516)
  * Dinh et al. (2017) Dinh L., Sohl-Dickstein J., Bengio S., 2017, Density estimation using Real NVP (arXiv:1605.08803)
  * Dome et al. (2022) Dome T., Fialkov A., Mocz P., Schäfer B. M., Boylan-Kolchin M., Vogelsberger M., 2022, Monthly Notices of the Royal Astronomical Society
  * Dome et al. (2023) Dome T., Fialkov A., Sartorio N., Mocz P., 2023, MNRAS, 525, 348
  * Dunstan et al. (2011) Dunstan R. M., Abazajian K. N., Polisensky E., Ricotti M., 2011, arXiv e-prints, p. arXiv:1109.6291
  * Font-Ribera et al. (2012) Font-Ribera A., et al., 2012, J. Cosmology Astropart. Phys., 2012, 059
  * Friedman & Hassan (2022) Friedman R., Hassan S., 2022, HIGlow: Conditional Normalizing Flows for High-Fidelity HI Map Modeling (arXiv:2211.12724)
  * Furlanetto et al. (2006) Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Phys. Rep., 433, 181
  * Gao et al. (2015) Gao L., Theuns T., Springel V., 2015, MNRAS, 450, 45
  * Ghara et al. (2017) Ghara R., Choudhury T. R., Datta K. K., Choudhuri S., 2017, MNRAS, 464, 2234
  * Giri & Schneider (2022) Giri S. K., Schneider A., 2022, Phys. Rev. D, 105, 083011
  * Hand et al. (2018) Hand N., Feng Y., Beutler F., Li Y., Modi C., Seljak U., Slepian Z., 2018, AJ, 156, 160
  * Hassan et al. (2019) Hassan S., Liu A., Kohn S., La Plante P., 2019, MNRAS, 483, 2524
  * Hassan et al. (2020) Hassan S., Andrianomena S., Doughty C., 2020, MNRAS, 494, 5761
  * Hassan et al. (2022) Hassan S., et al., 2022, ApJ, 937, 83
  * Hložek et al. (2015) Hložek R., Grin D., et al., 2015, PRD, 91, 103512
  * Ho et al. (2021) Ho M.-F., Bird S., Garnett R., 2021, MNRAS, 507, 704
  * Hu et al. (2000) Hu W., Barkana R., et al., 2000, Phys. Rev. Lett., 85
  * Iršič et al. (2023) Iršič V., et al., 2023, arXiv e-prints, p. arXiv:2309.04533
  * Jones et al. (2018) Jones M. G., Haynes M. P., Giovanelli R., Moorman C., 2018, MNRAS, 477, 2
  * Jones et al. (2021) Jones D., Palatnick S., Chen R., Beane A., Lidz A., 2021, ApJ, 913, 7
  * Kaiser (1987) Kaiser N., 1987, MNRAS, 227, 1
  * Kakiichi et al. (2017) Kakiichi K., et al., 2017, MNRAS, 471, 1936
  * Kim & Marsh (2016) Kim J. E., Marsh D. J. E., 2016, Phys. Rev. D, 93, 025027
  * Kingma & Dhariwal (2018) Kingma D. P., Dhariwal P., 2018, Glow: Generative Flow with Invertible 1x1 Convolutions (arXiv:1807.03039)
  * Lagos et al. (2014) Lagos C. D. P., Baugh C. M., Zwaan M. A., Lacey C. G., Gonzalez-Perez V., Power C., Swinbank A. M., van Kampen E., 2014, MNRAS, 440, 920
  * Lidz & Hui (2018) Lidz A., Hui L., 2018, Phys. Rev. D, 98, 023011
  * Lin et al. (2020) Lin X., Cai Z., Li Y., Krolewski A., Ferraro S., 2020, ApJ, 905, 176
  * List et al. (2020) List F., Elahi P. J., Lewis G. F., 2020, ApJ, 904, 153
  * Liu & Shaw (2020) Liu A., Shaw J. R., 2020, PASP, 132, 062001
  * Liu & Slatyer (2018) Liu H., Slatyer T. R., 2018, arXiv e-prints, p. arXiv:1803.09739
  * Madelung (1927) Madelung E., 1927, Zeitschrift fur Physik, 40, 322
  * Marsh (2016) Marsh D. J. E., 2016, Physics Reports, 643, 1
  * May & Springel (2021) May S., Springel V., 2021, MNRAS, 506, 2603
  * McKee & Krumholz (2009) McKee C. F., Krumholz M. R., 2009, The Astrophysical Journal, 709, 308
  * Mocz et al. (2018) Mocz P., Lancaster L., Fialkov A., Becerra F., Chavanis P.-H., 2018, PRD, 97, 083519
  * Mocz et al. (2019) Mocz P., et al., 2019, PRL, 123, 141301
  * Mocz et al. (2023) Mocz P., et al., 2023, MNRAS, 521, 2608
  * Muñoz & Loeb (2018) Muñoz J. B., Loeb A., 2018, arXiv e-prints, p. arXiv:1802.10094
  * Nebrin et al. (2019) Nebrin O., Ghara R., Mellema G., 2019, J. Cosmology Astropart. Phys., 2019, 051
  * Newburgh et al. (2016) Newburgh L. B., et al., 2016, in Hall H. J., Gilmozzi R., Marshall H. K., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 9906, Ground-based and Airborne Telescopes VI. p. 99065X (arXiv:1607.02059), doi:10.1117/12.2234286
  * Ni et al. (2019) Ni Y., Wang M.-Y., Feng Y., Di Matteo T., 2019, MNRAS, 488, 5551
  * Noterdaeme et al. (2009) Noterdaeme P., Petitjean P., Ledoux C., Srianand R., 2009, A&A, 505, 1087
  * Noterdaeme et al. (2012) Noterdaeme P., et al., 2012, A&A, 547, L1
  * Oman (2022) Oman K. A., 2022, MNRAS, 509, 3268
  * Paduroiu (2022) Paduroiu S., 2022, Universe, 8
  * Papamakarios et al. (2021) Papamakarios G., Nalisnick E., Rezende D. J., Mohamed S., Lakshminarayanan B., 2021, Journal of Machine Learning Research, 22, 1
  * Peccei & Quinn (1977) Peccei R. D., Quinn H. R., 1977, Phys. Rev. Lett., 38, 1440
  * Pérez-Ràfols et al. (2018) Pérez-Ràfols I., et al., 2018, MNRAS, 473, 3019
  * Péroux et al. (2005) Péroux C., Dessauges-Zavadsky M., D’Odorico S., Sun Kim T., McMahon R. G., 2005, MNRAS, 363, 479
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Pober et al. (2013) Pober J. C., et al., 2013, AJ, 145, 65
  * Pober et al. (2014) Pober J. C., et al., 2014, ApJ, 782, 66
  * Pritchard & Loeb (2012) Pritchard J. R., Loeb A., 2012, Reports on Progress in Physics, 75, 086901
  * Puchwein et al. (2023) Puchwein E., et al., 2023, MNRAS, 519, 6162
  * Rahmati et al. (2015) Rahmati A., Schaye J., Bower R. G., Crain R. A., Furlong M., Schaller M., Theuns T., 2015, MNRAS, 452, 2034
  * Sadler et al. (2020) Sadler E. M., et al., 2020, MNRAS, 499, 4293
  * Santos et al. (2015) Santos M., et al., 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14). p. 19 (arXiv:1501.03989), doi:10.22323/1.215.0019
  * Sarkar et al. (2016) Sarkar D., Bharadwaj S., Anathpindika S., 2016, MNRAS, 460, 4310
  * Schive et al. (2014) Schive H.-Y., Chiueh T., Broadhurst T., 2014, Nature Physics, 10, 496
  * Schneider (2018) Schneider A., 2018, Phys. Rev. D, 98, 063021
  * Schneider et al. (2012) Schneider A., Smith R. E., Macciò A. V., Moore B., 2012, MNRAS, 424, 684
  * Sheth & Tormen (2002) Sheth R. K., Tormen G., 2002, MNRAS, 329, 61
  * Sitwell et al. (2014) Sitwell M., Mesinger A., Ma Y.-Z., Sigurdson K., 2014, MNRAS, 438, 2664
  * Songaila & Cowie (2010) Songaila A., Cowie L. L., 2010, ApJ, 721, 1448
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
  * Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
  * Switzer et al. (2013) Switzer E. R., et al., 2013, MNRAS, 434, L46
  * Tegmark et al. (2004) Tegmark M., et al., 2004, ApJ, 606, 702
  * Theuns (2021) Theuns T., 2021, MNRAS, 500, 2741
  * Tinker et al. (2010) Tinker J. L., Robertson B. E., Kravtsov A. V., Klypin A., Warren M. S., Yepes G., Gottlöber S., 2010, ApJ, 724, 878
  * Villaescusa-Navarro (2018) Villaescusa-Navarro F., 2018, Pylians: Python libraries for the analysis of numerical simulations, Astrophysics Source Code Library, record ascl:1811.008 (ascl:1811.008)
  * Villaescusa-Navarro et al. (2014) Villaescusa-Navarro F., Viel M., Datta K. K., Choudhury T. R., 2014, J. Cosmology Astropart. Phys., 2014, 050
  * Villaescusa-Navarro et al. (2015) Villaescusa-Navarro F., Bull P., Viel M., 2015, ApJ, 814, 146
  * Villaescusa-Navarro et al. (2018) Villaescusa-Navarro F., et al., 2018, ApJ, 866, 135
  * Villaescusa-Navarro et al. (2021) Villaescusa-Navarro F., et al., 2021, ApJ, 915, 71
  * Vogelsberger et al. (2013) Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L., 2013, MNRAS, 436, 3031
  * Wang et al. (2021) Wang Z., et al., 2021, ApJ, 907, 4
  * Weinberg et al. (2015) Weinberg D. H., Bullock J. S., et al., 2015, Proceedings of the National Academy of Sciences, 112
  * Weinberger et al. (2020) Weinberger R., Springel V., Pakmor R., 2020, ApJS, 248, 32
  * Wyithe & Loeb (2008) Wyithe J. S. B., Loeb A., 2008, MNRAS, 383, 606
  * Yang et al. (2020) Yang J., et al., 2020, ApJ, 904, 26
  * Zafar et al. (2013) Zafar T., Péroux C., Popping A., Milliard B., Deharveng J. M., Frank S., 2013, A&A, 556, A141
  * Zamudio-Fernandez et al. (2019) Zamudio-Fernandez J., Okan A., Villaescusa-Navarro F., Bilaloglu S., Derin Cengiz A., He S., Perreault Levasseur L., Ho S., 2019, arXiv e-prints, p. arXiv:1904.12846