  The joint likelihood of observable cluster signals reflects the astrophysical
evolution of the coupled baryonic and dark matter components in massive halos,
and its knowledge will enhance cosmological parameter constraints in the coming
era of large, multi-wavelength cluster surveys. We present a computational
study of intrinsic covariance in cluster properties using halo populations
derived from Millennium Gas Simulations (MGS). The MGS are re-simulations of
the original 500 Mpc/h Millennium Simulation performed with gas dynamics under
two different physical treatments: shock heating driven by gravity only (GO)
and a second treatment with cooling and preheating (PH). We examine
relationships among structural properties and observable X-ray and
Sunyaev-Zel'dovich (SZ) signals for samples of thousands of halos with M_200 >
5 \times 10^{13} Msun/h and z < 2. While the X-ray scaling behavior of PH model
halos at low-redshift offers a good match to local clusters, the model exhibits
non-standard features testable with larger surveys, including weakly running
slopes in hot gas observable--mass relations and ~10% departures from
self-similar redshift evolution for 10^14 Msun/h halos at redshift z ~ 1. We
find that the form of the joint likelihood of signal pairs is generally
well-described by a multivariate, log-normal distribution, especially in the PH
case which exhibits less halo substructure than the GO model. At fixed mass and
epoch, joint deviations of signal pairs display mainly positive correlations,
especially the thermal SZ effect paired with either hot gas fraction
(r=0.88/0.69 for PH/GO at z=0) or X-ray temperature (r=0.62/0.83). We discuss
halo mass selection by signal pairs, and find a minimum mass scatter of 4% in
the \PH model by combining thermal SZ and gas fraction measurements.
