  Intensity mapping (IM) of spectral lines has the potential to revolutionize
cosmology by increasing the total number of observed modes by several orders of
magnitude compared to the cosmic microwave background (CMB) anisotropies. In
this paper, we consider IM of neutral hydrogen (HI) in the redshift range $0
\lesssim z \lesssim 3$ employing a halo model approach where HI is assumed to
follow the distribution of dark matter (DM) halos. If a portion of the DM is
composed of ultralight axions then the abundance of halos is changed compared
to cold dark matter below the axion Jeans mass. With fixed total HI density,
$\Omega_{\rm HI}$, assumed to reside entirely in halos, this effect introduces
a scale-independent increase in the HI power spectrum on scales above the axion
Jeans scale, which our model predicts consistent with N-body simulations.
Lighter axions introduce a scale-dependent feature even on linear scales due to
its suppression of the matter power spectrum near the Jeans scale. We use the
Fisher matrix formalism to forecast the ability of future HI surveys to
constrain the axion fraction of DM and marginalize over astrophysical and model
uncertainties. We find that a HIRAX-like survey is a very reliable IM survey
configuration, being affected minimally by uncertainties due to non-linear
scales, while the SKA1MID configuration is the most constraining as it is
sensitive to non-linear scales. Including non-linear scales and combining a
SKA1MID-like IM survey with the Simons Observatory CMB, the benchmark "fuzzy
DM" model with $m_a = 10^{-22}\text{ eV}$ can be constrained at the 10% level.
For lighter ULAs this limit improves below 1%, and allows the possibility to
test the connection between axion models and the grand unification scale across
a wide range of masses.
