  We use an analytical forward model based on perturbation theory to predict
the neutral hydrogen (HI) overdensity maps at low redshifts. We investigate its
performance by comparing it directly at the field level to the simulated HI
from the IllustrisTNG simulation TNG300-1 ($L=205\ h^{-1}$ Mpc), in both real
and redshift space. We demonstrate that HI is a biased tracer of the underlying
matter field and find that the cubic bias model describes the simulated HI
power spectrum to within 1% up to $k=0.4 \;(0.3) \,h\,{\rm Mpc}^{-1}$ in real
(redshift) space at redshifts $z=0,1$. Looking at counts in cells, we find an
excellent agreement between the theory and simulations for cells as small as 5
$h^{-1}$ Mpc. These results are in line with expectations from perturbation
theory and they imply that a perturbative description of the HI field is
sufficiently accurate given the characteristics of upcoming 21cm intensity
mapping surveys. Additionally, we study the statistical properties of the model
error - the difference between the truth and the model. We show that on large
scales this error is nearly Gaussian and that it has a flat power spectrum,
with amplitude significantly lower than the standard noise inferred from the HI
power spectrum. We explain the origin of this discrepancy, discuss its
implications for the HI power spectrum Fisher matrix forecasts and argue that
it motivates the HI field-level cosmological inference. On small scales in
redshift space we use the difference between the model and the truth as a proxy
for the Fingers-of-God effect. This allows us to estimate the nonlinear
velocity dispersion of HI and show that it is smaller than for the typical
spectroscopic galaxy samples at the same redshift. Finally, we provide a simple
prescription based on the perturbative forward model which can be used to
efficiently generate accurate HI mock data, in real and redshift space.
