  Line intensity mapping (LIM) is an emerging technique with a unique potential
to probe a wide range of scales and redshifts. Realizing the full potential of
LIM, however, relies on accurate modeling of the signal. We introduce an
extended halo model for the power spectrum of intensity fluctuations of CO
rotational lines and [CII] fine transition line in real space, modeling
nonlinearities in matter fluctuations and biasing relation between the line
intensity fluctuations and the underlying dark matter distribution. We also
compute the stochastic contributions beyond the Poisson approximation using the
halo model framework. To establish the accuracy of the model, we create the
first cosmological-scale simulations of CO and [CII] intensity maps,
\textsf{MithraLIMSims}, at redshifts $0.5 \leq z\leq6$, using halo catalogs
from Hidden-Valley simulations, and painting halos according to
mass-redshift-luminosity relations for each line. We show that at $z=1$ on
scales $k_{\rm max} \lesssim 0.8 \ {\rm Mpc}^{-1}h$, the model predictions of
clustering power (with only two free parameters) are in agreement with the
measured power spectrum at better than 5\%. At higher redshift of $z=4.5$, this
remarkable agreement extends to smaller scale of $ k_{\rm max} \lesssim 2 \
{\rm Mpc}^{-1}h$. Furthermore, we show that on large scales, the stochastic
contributions to CO and CII power spectra are non-Poissonian, with amplitudes
reproduced reasonably well by the halo model prescription. Lastly, we assess
the performance of the theoretical model of the baryon acoustic oscillations
(BAO) and show that hypothetical LIM surveys probing CO lines at $z=1$, that
can be deployed within this decade, will be able to make a high significance
measurement of the BAO. On a longer time scale, a space-based mission probing
[CII] line can uniquely measure the BAO on a wide range of redshifts at an
unprecedented precision.
