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 magneto-hydrodynamical 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 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 21 cm 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.
