Not yet reviewed by Pith; the record is open.

The tightest cosmological constraints currently available are obtained by combining complementary data sets. When combining correlated data sets, various astrophysical biases that affect the measurements must be identified and treated. There are numerous such biases, and they are often intricately related with one another via complex astrophysical effects, making them difficult to characterize analytically. Consequently, a simulation with multiple components implemented coherently is required to investigate these biases simultaneously and as a whole. In this work, a suite of simulated extragalactic skies is presented, including maps and/or catalogues of cosmic microwave background (CMB) lensing, thermal and kinetic Sunyaev-Zel'dovich (tSZ/kSZ) effects, cosmic infrared background (CIB), radio sources, galaxy overdensity and galaxy weak lensing. Each of these probes is implemented in the lightcone using halo catalogues and/or particles from the Multidark-Planck2 (MDPL2) N-body simulation, and the modelling is calibrated using hydrodynamic simulations and publicly available data. The auto- and cross-spectra of the individual probes, as well as the cross-spectra between the observables, are shown to be consistent with theoretical models and measurements from data. The simulation is shown to have a wide range of applications, including forecasting, pipeline testing, and evaluating astrophysical biases in cross-correlation studies. It is further demonstrated that the simulation products produced in this work have sufficient accuracy to recover the input cosmology when subjected to a full cosmological analysis and are ready for application in real-world analyses for ongoing and future surveys.