  Intensity mapping -- the large-scale mapping of selected spectral lines
without resolving individual sources -- is quickly emerging as an efficient way
to conduct large cosmological surveys. Multiple surveys covering a variety of
lines (such as the hydrogen 21cm hyperfine line, CO rotational lines, and [CII]
fine structure lines, among others) are either observing or will soon be
online, promising a panchromatic view of our Universe over a broad redshift
range. With multiple lines potentially covering the same volume,
cross-correlations have become an attractive prospect, both for probing the
underlying astrophysics and for mitigating observational systematics. For
example, cross correlating 21cm and [CII] intensity maps during reionization
could reveal the characteristic scale of ionized bubbles around the first
galaxies, while simultaneously providing a convenient way to reduce independent
foreground contaminants between the two surveys. However, many of the desirable
properties of cross-correlations in principle emerge only under ideal
conditions, such as infinite ensemble averages. In this paper, we construct an
end-to-end pipeline for analyzing intensity mapping cross-correlations,
enabling instrumental effects, foreground residuals, and analysis choices to be
propagated through Monte Carlo simulations to a set of rigorous error
properties, including error covariances, window functions, and full probability
distributions for power spectrum estimates. We use this framework to critically
examine the applicability of simplifying assumptions such as the independence
and Gaussianity of power spectrum errors. As worked examples, we forecast the
sensitivity of near-term and futuristic 21cm-[CII] cross-correlation
measurements, providing recommendations for survey design.
