  Intensity mapping surveys will provide access to a coarse view of the cosmic
large-scale structure in unprecedented large volumes at high redshifts. Given
the large fractions of the sky that can be efficiently scanned using emission
from cosmic neutral hydrogen (HI), intensity mapping is ideally suited to probe
a wide range of density environments and hence to constrain cosmology and
fundamental physics. To efficiently extract information from 21cm intensities
beyond average, one needs non-Gaussian statistics that capture large deviations
from mean HI density. Counts-in-cells statistics are ideally suited for this
purpose, as the statistics of matter densities in spheres can be predicted
accurately on scales where their variance is below unity.
  We use a large state-of-the-art magneto-hydrodynamic simulation from the
IllustrisTNG project to determine the relation between neutral hydrogen and
matter densities in cells. We demonstrate how our theoretical knowledge about
the matter PDF for a given cosmology can be used to extract a
parametrisation-independent HI bias function from a measured neutral hydrogen
PDF. When combining the predicted matter PDFs with a simple bias fit to the
simulation, we obtain a prediction for neutral hydrogen PDFs at a few percent
accuracy at scale R=5 Mpc/h from redshift z=5 to z=1. Furthermore, we find a
density-dependent HI clustering signal that is consistent with theoretical
expectations and could allow for joint constraints of HI bias and the amplitude
of matter fluctuations or the growth of structure.
