  The 21 cm intensity mapping experiments promise to obtain the large-scale
distribution of HI gas at the post-reionization epoch. In order to reveal the
underlying matter density fluctuations from the HI mapping, it is important to
understand how HI gas traces the matter density distribution. Both nonlinear
halo clustering and nonlinear effects modulating HI gas in halos may determine
the scale below which the HI bias deviates from linearity. We employ three
approaches to generate the mock HI density from a large-scale N-body simulation
at low redshifts, and demonstrate that the assumption of HI linearity is valid
at the scale corresponding to the first peak of baryon acoustic oscillations,
but breaks down at $k \gtrsim 0.1\,h\, {\rm Mpc}^{-1}$. The nonlinear effects
of halo clustering and HI content modulation counteract each other at small
scales, and their competition results in a model-dependent "sweet-spot"
redshift near $z$=1 where the HI bias is scale-independent down to small
scales. We also find that the linear HI bias scales approximately linearly with
redshift for $z\le 3$.
