{'fraction_non_alphanumeric': 0.03428564388759532, 'fraction_numerical': 0.018454866885956808, 'mean_word_length': 4.616273438040545, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Websky Extragalactic CMB Simulations

We present a new pipeline for the efficient generation of synthetic observations of the extragalactic microwave sky, tailored to large ground-based CMB experiments such as the Simons Observatory, Advanced ACTPol, SPT-3G, and CMB-S4. Such simulated observations are a key technical challenge in cosmology because of the dynamic range and accuracy required. The first part of the pipeline generates a random cosmological realization in the form of a dark matter halo catalog and matter displacement field, as seen from a given position. The halo catalog and displacement field are modeled with ellipsoidal collapse dynamics and Lagrangian perturbation theory, respectively. In the second part, the cosmological realization is converted into a set of intensity maps over the range 10 - 10^3 GHz using models based on existing observations and hydrodynamical simulations. These maps include infrared emission from dusty star forming galaxies (CIB), Comptonization of CMB photons by hot gas in groups and clusters through the thermal Sunyaev-Zel'dovich effect (tSZ), Doppler boosting by Thomson scattering of the CMB by bulk flows through the kinetic Sunyaev-Zel'dovich effect (kSZ), and weak gravitational lensing of primary CMB anisotropies by the large-scale distribution of matter in the universe. After describing the pipeline and its implementation, we present the Websky maps, created from a realization of the cosmic web on our past light cone in the redshift interval 0<z<4.6 over the full-sky and a volume of ~(600 Gpc/h)^3 resolved with ~10^12 resolution elements. The Websky maps and halo catalog are publicly available at mocks.cita.utoronto.ca/websky.

Previous simulations modeled the different components with varying levels of accuracy, depending on the application, such as extragalactic point sources [28][29][30], the tSZ effect in small regions (from individual clusters subtending a few arcmin to fields of view of ∼ 100 sq deg.) [31][32][33][34][35][36][37][38][39][40][41][42][43], the CIB [44,45], the kSZ [46], and CMB lensing [47][48][49]. State-of-the-art massively parallel cosmological simulations with hydrodynamics and subgrid feedback prescriptions still require considerable computational resources to generate even a single realization, and are currently limited to a maximum size of ∼ 1 Gpc 3 (although generally much smaller) [50][51][52][53][54][55][56][57]. Full-sky simulations that map out a volume corresponding to our past light cone therefore require approximate large-scale structure realizations without explicit hydrodynamic simulation of individual objects. Rather, small scale observations and hydrodynamic simulations can be used to calibrate efficient large-scale models by a halo model [58] or the dark matter-only distribution [59]. Large-scale cosmological N-body simulations focused on halo catalogs and the overall cosmic web of dark matter can be done much more quickly than high-resolution galaxy formation simulations, and were used as the input for a model to simulate the extragalactic microwave sky in [60] and [61], including synchrotron and dust in galaxies, and lensing, with the resulting maps being used extensively in modeling of data from ACT and other CMB experiments to the present day. This approach of applying hydrodynamical results to N-body simulations was further validated in [58], where it was demonstrated that explicitly replacing pressure profiles in hydrodynamic simulations with the spherically-averaged mean for halos of a given mass and redshift produces similar results.
An efficient new pipeline we have developed that builds on this previous work, and the publicly available sky maps generated with it, are the focus of this paper. The new pipeline combines (1) a parallel on-the-fly lightcone algorithm for fast halo catalogs with (2) a halo-field model for the observables, to generate simulations of the extragalactic sky from microwave to far infrared bands much more efficiently than previously possible. Among the available accelerated methods for generating realizations of large scale structure [62][63][64][65][66][67], we have chosen to use the mass-Peak Patch approach [68,69] for this work. In §2 we describe the first step in the pipeline, efficient generation of a realization of LSS on the past light cone. In §3 we describe the second step, going from the LSS realization to maps of the sky at different wavelengths, from microwave to the far infrared. The results of the simulations compared to recent observational data are given in §4, including angular auto and cross power spectra of the CIB, tSZ, kSZ, and lensing signals. We conclude with a summary and discussion of future directions in §5. The following flat ΛCDM cosmological parameters, consistent with Planck 2018 [70], are used throughout this paper: (Ω m , Ω b , σ 8 , n s , h, τ) = (0.31, 0.049, 0.81, 0.965, 0.68, 0.055).

Large Scale Structure Simulations
The sky fraction, redshift coverage, and mass resolution required for accurate extragalactic mocks of current and next-generation CMB observations demand exceedingly large simulations. For example, a single-box full-sky lightcone simulation out to z = 5, capable of resolving dark matter halos above a minimum mass of ∼1× 10 12 M /h with 20 resolution elements each, requires at least ∼12,000 3 total resolution elements in a box size of ∼10 Gpc/h. This will remain computational infeasible for hydrodynamical methods over the next decade, and a single realization would be currently achievable through full N-body simulation only at the expense of many millions of CPU hours -a considerable investment for all but the largest computational grants. One approach to alleviate this computational bottleneck is to instead simulate a smaller volume with either hydrodynamical or N-body methods and replicate it as needed to fill the sky fraction and redshift requirements. Unfortunately, this repeats structures along the line of sight, does not reproduce the statistics of rare objects, and fails to capture the effects of the large-scale density and velocity fluctuations that are important for accurate mocks of gravitational lensing and the kSZ effect.
The demand for rapid generation of these type of simulations has therefore resulted in the development of many 'approximate' methods for generating realizations of large scale structure [62][63][64][65][66][67][68][71][72][73][74]. Many of these simulation methods were recently featured and discussed in a comparison project within the Euclid collaboration focused on determining their viability for estimating covariance matrices [75][76][77]. This work provided an extensive analysis of the halo correlation function, halo power spectrum, and halo bispectrum, and their individual covariances, from a set of 300 approximate simulations at z = 1, in comparison to the results of full N-body. To generate the dark matter halo catalogues for this work we use the mass-Peak Patch approach, recently described and validated in detail in [68] and included in the Euclid comparison project. The mass-Peak Patch approach has been extensively compared to more expensive N-body simulation runs and alternative approximate methods, with highly-satisfactory results at the simulation resolutions required for this work. Individual peak-to-group comparisons show good agreement for high-mass, tightly bound groups, with growing scatter for lower masses and looser binding. The final state (Eulerian) spatial distribution of peak patches and N-body clusters, and halo velocities, have been shown to be satisfyingly close. As well, higher order halo statistics such as the halo bispectrum have been shown to be accurately reproduced when considering abundance-matched samples [77] such as presented here. We use mass-Peak Patch in this work because it is computationally inexpensive, generates light cone catalogs on the fly with no explicit time-stepping, and the method fared well in the Euclid comparison project and subsequent validations. Using the same sky model applied in this work on a large scale structure realization generated using a different approximate simulation method should produce statistically similar extragalactic maps.
The cosmological realizations presented here include halos and the material exterior to halos, the latter of which we refer to as the 'field' component. The halo properties are determined using the mass-Peak Patch approach [68,69], while the dynamics of the field component, mainly of applicability for creating kSZ and lensing maps, is based on second order Lagrangian perturbation theory (2LPT) [78]. We provide a brief summary of the mass-Peak Patch method below, but refer the reader to [68] for detailed explanations and validation results. The light cone pipeline can be separated into four main subprocesses: 1. Generation of a linear random Gaussian density field, δ , and the corresponding 2LPT displacement vectors s (1) and s (2)

Initial Conditions & Domain Decomposition
The initial conditions for a mass-Peak Patch simulation, as is standard for the vast majority of largescale structure simulations, are in the form of a cubic periodic lattice of any specified physical size L and one dimensional resolution n. The total number of lattice sites, or equivalently particles, in the simulation is then N = n 3 , initially uniformly distributed in a total volume of V = L 3 . The uniform periodic lattice allows many operations to be performed in the Fourier domain to accelerate computation of convolutions and derivatives, where here we define the relation of a Lagrangian field f (q) and its Fourier transform f (k) as f (q) = (2π) −3 dke ik·q f (k).
To generate the initial conditions on a periodic lattice, random perturbations are realized by first generating a white noise field by drawing each lattice value from a Gaussian distribution with a unit variance and a mean of zero. It is then convolved with the linear matter power spectrum to obtain the linear density contrast δ (k). The first and second order linear displacements, describing the displacement of matter as a function of time and Lagrangian position, are obtained using 2LPT, requiring the storage of seven values at each lattice site: {δ (q), s (1) (q), s (2) (q)}. All quantities are expressed in terms of their values at z = 0, and scaled back in time by the appropriate linear growth factor where required, unless otherwise noted.
The parallelization of the subsequent peak finding algorithm and homogeneous ellipsoid calculations are based on a domain decomposition consisting of overlapping cubic 'tiles', which greatly reduces the need for communication and collective operations in the peak finding and exclusion algorithms. This choice allows maximum flexibility in separating the cosmological fluctuations into long and short-range components, with the possibility for a coarser resolution for long wavelength modes and adaptive resolution of smaller structures. For the full-sky simulations presented in this paper, uniform spatial resolution over a spherical region is required, and so the initial conditions are realized at uniform resolution on the periodic cubic lattice. Convolutions that require collective communication, such as calculation of displacements from the linear density field using perturbation theory, are done globally in parallel at the native resolution of the realization with the use of the slab-decomposed FFTW3 library 3 . Depending on memory constraints, the density and displacement fields are stored either in main memory with a slab-decomposition, or as segmented files on disk for later processing.
Once initial condition generation is complete, exchange is performed between processes to go from a slab to cubic domain decomposition with N 3 t tiles. For example, an eight-process run would start by realizing the linear density contrast and displacements in eight slabs, one for each processor. After this initial step, the fields would be rearranged on the processes such that there are eight cubic tiles in a 2 x 2 x 2 configuration, N t = 2. Because the subsequent peak finding in each tile is done in parallel without any communication between tiles, including smoothing the density and displacements with a spherically-symmetric filter, it is necessary to retain information about the field in the nearby adjacent tiles through the use of a buffer region. The buffer region is taken to have a width given by the largest possible halo radius in Lagrangian coordinates. For example, in a Gpc-size region at z = 0, one can expect to find halos as large as several times 10 15 M , corresponding to a Lagrangian radius of about 40 Mpc. The smoothing can then be done locally using a threaded FFT on the tile with no further domain decomposition. For further details of the parallelization scheme and computational requirements see Appendix A of [68].

Ellipsoidal Collapse
The halo catalogs are based on the parallel implementation of the mass-Peak Patch approach [68], wherein ellipsoidal collapse is used to determine locations and masses of potential collapsed objects (halos). The mass-Peak Patch approach is a Lagrangian space halo finder that associates halos with the largest regions that have just collapsed by a given redshift or, equivalently, distance from the observer. The determination of whether any given region will have collapsed or not is made by approximating it as a homogeneous ellipsoid, the fate of which is determined completely by the principal axes of the deformation tensor of the linear displacement field (i.e. the strain) averaged over the region. In principle, the process of finding these local mass peaks would involve measuring the strain at every point in space, smoothed on every scale. However, experimentation has shown that equivalent results can be obtained by measuring the strain around density peaks found on a range of scales. This is not to say that a halo found on a given scale corresponds to a peak in the density smoothed on that scale, however, which is only the case when the strain is isotropic and the collapse is spherical. Thus, the use of density peaks as centers for strain measurements and ellipsoidal collapse calculations in the algorithm is only an optimization, to avoid wasting computations measuring the properties of regions of Lagrangian space that will not collapse in the first place.
The first step in identifying the possible site of a halo is therefore to identify peaks in the linear density field. This is done by smoothing the field on a series of top hat filters of different scales, which we will refer to as the 'filter bank'. The filter bank can consist of logarithmically-spaced filter radii R, or linear of logarithmic spacing in σ (R), with optimal filter spacings to maximize both accuracy and efficiency presented in [68], from a minimum radius of R f,min = 2a latt , where a latt is the lattice spacing, to a maximum radius of the largest halo expected to be found in the simulated volume (R f,max = 36 Mpc at z = 0). Starting with the largest smoothing scale, a local FFT is performed on the local tile in order to convolve the field with the top hat smoothing kernel W (kR) = 3 j 1 (kR)/(kR) where j 1 is the spherical Bessel function of order one. All points that are density peaks on that filter scale with a density contrast above a threshold determined from spherical collapse, roughly 1.5, are added to the list of candidate halo centers with a 'peak radius' given by the scale at the filter scale, R pk = R f . This process is repeated for each smoothing scale in the filter bank until a full list of candidate halo centers has been compiled.
Once the candidate list is completed, measurements are made by averaging the linear strain, e i j ≡ −(∂ s (1) i /∂ q j + ∂ s (1) j /∂ q i )/2, over spheres with Lagrangian radii R, centered on each halo candidate. The measurements are started from some initial radius chosen to be a fixed multiple of the smoothing scale on which the candidate center was found to be a density peak R init = f init R pk . If the region averaged over the initial radius is determined to collapse, the measurements are performed at progressively larger radii until the region no longer collapses, while if the region averaged over the initial radius does not collapse, the measurements are performed at decreasing radii until the region collapses, or is equal to a radius smaller than a lattice size. The mean strain within a given measurement radius is diagonalized to obtain its eigenvalues, λ i = −(δ (1) /3)(1 + c i ), and eigenvectors, n i , with e i j = ∑ k λ kn i kn j k , c 1 = p + 3e, c 2 = −2p, and c 3 = p − 3e, where the linear density contrast and strain, δ (1) and e i j , respectively, are averaged over the smoothing radius R. The redshift at which the corresponding homogeneous ellipsoid would collapse is calculated from the following system of ordinary differential equations:ẍ where x i (t) = R i (t)/R m are the scale factors in each of the three principal axes of the ellipsoid, a is the background scale factor, δ = a 3 /(x 1 x 2 x 3 ) − 1 is the evolving non-linear density contrast of the homogeneous ellipsoid, δ lin = δ (1) D(a), and The parameters b i are defined by the Newtonian potential of an isolated homogeneous ellipsoid in a coordinate system aligned with its principal axes: Φ(x) = 4πGρ ∑ i b i x 2 i /2, and depend only on its evolving shape. Therefore, the first term in the brackets on the right hand side of equation (2.1) is proportional to the gravitational acceleration from material inside of the ellipsoid, and is exact, while the second term represents that from external tidal fields, and is approximated by the linear solution, which only depends on the mean linear strain within R, through the c i [69]. The parameters b i and c i can be considered to encode the anisotropy of the internal and external forces, respectively. Since the shape of the ellipsoid is time dependent, while the shape of the linear tidal field is not, b i is time dependent while c i is not. Note that in the spherical case, x 1 = x 2 = x 3 = x implies b i = 1, c i = 0, andẍ/x =ä/a − Ω m H 2 δ /2, which is the usual equation of motion for spherical collapse. The initial conditions at time t 0 for these equations are obtained using perturbation theory. Each axis of the ellipsiod is evolved until it reaches a critical radius x eq,i = f r,i a during its collapse, after which it is 'frozen in' at that value. Since the eigenvalues of the ellipsoid satisfy λ v3 ≥ λ v2 ≥ λ v1 , the 1-axis will be last to collapse, and it is at this point when a peak is considered virialized. The radial freezeout factors f r,3 and f r,2 were chosen to be 0.171, as this corresponds to the standard top hat virial density contrast of 200. The final axis was frozen out at f r,1 = 0.01, indicating complete collapse along the final axis.
For each candidate peak we find the largest radius for which a homogeneous ellipsoid with the measured mean strain in the sphere contained within the radius would collapse by the redshift of interest, z = z coll (δ m , e m , p m ), with the collapse redshift function z coll stored in a precomputed table for computational efficiency. If a candidate peak has no radius for which a homogeneous ellipsoid with the measured strain would have collapsed, then that point is discarded. Each candidate position and radius is then stored as a peak patch. We then proceed down through the filter bank and repeat this procedure for each scale, resulting in a list of peak patches which we refer to as the 'un-merged catalog'.

Exclusion and Displacements
Exclusion is essential to avoid double counting of matter in halos, since distinct halos do not overlap, by definition. Because neighboring and overlapping regions in Lagrangian space undergo ellipsoidal collapse in our calculations, we must use some practical algorithm for ensuring that exclusion is taken into account. The approach here is essentially unchanged from "binary exclusion" as defined in section 3.2 of [69].
Binary exclusion starts from a list of candidate peak patches sorted by mass or, equivalently, Lagrangian peak patch radius. For each patch we consider every other less massive patch that overlaps it. If the center of the smaller patch is inside the large one, then that patch is removed from the list. If the center of the smaller patch is outside of the larger one, then the volume of overlap of the two intersecting spheres is calculated from the standard equations describing the volume of their overlapping spherical caps. Considering the plane perpendicular to the separation vector between them, we subtract any volume beyond this plane from each of the halos. The overlapping mass from all neighbouring halos is accounted for before reducing their volumes. This process is repeated until the least massive remaining patch in the list is reached.
We use Lagrangian perturbation theory (LPT) to move halos and the matter outside them from their initial Lagrangian positions to their final Eulerian positions on the light cone. The first step is to compute the 1LPT and 2LPT displacements for each of the mass elements in the simulation, which is done globally in parallel when generating the initial conditions. The 1LPT displacements are obtained from s (1) (k) = −iδ (1) (k)k/k 2 and the second order displacements are obtained with s (2) (k) = −iδ (2) (1) ii (q)e (1) j j (q) − e (1) i j (q)e (1) i j (q)]. In practice, we perform this sum by calculating the elements of the strain tensor in harmonic space, e (1) i j (k) = −k i k j δ (1) (k)/k 2 , and Fourier transforming to configuration space and performing the multiplications necessary sequentially, which maintains the memory usage of the code at seven floats per mass element (δ (1) , s (1) , s (2) ), at the expense of only one additional Fourier transform. During the adaptive measurement step to determine the largest collapsing region around candidate peak patch centers, we also keep track of the mean displacements, s (1) and s (2) , so that each unmerged object in the catalog is associated with its mean LPT displacement.

Light Cone
The algorithm we have presented here can be applied at fixed time or from a fixed vantage point -i.e., on the past light cone of the observer, which is the case for all the simulations we present here. Common choices for the observer are either at the center of the periodic simulation volume, with side λ box , or placed at a corner to consider only one octant of the sky. Every comoving cell in the simulation volume, at position q c , corresponds to a fixed redshift given a cosmological model, z c = z(χ c ), where χ c = |q c |. Each halo and field mass element is moved to its Eulerian position using , which is an excellent approximation [78].

Sky Model
The simulated extra-galactic microwave sky maps presented in this work are generated from the cosmological realizations by projecting, along the line of sight to each pixel and frequency, the emission from dusty star forming galaxies (CIB), Comptonization of CMB photons by hot gas in groups and clusters (tSZ), Doppler boosting by Thomson scattering of the CMB by bulk flows (kSZ), and weak lensing of the CMB by the intervening matter distribution. In §3.1 we first describe the general use of diffuse halo profiles, halo occupation distribution models, and treatment of field particles used, while §3.2 details the specific astrophysical models applied.

General Approach
The starting point for generating sky maps is the large scale structure along the observed lightcone, i.e. the mass and positions of halos and the LPT displacement field. For lensing and SZ secondaries, we use continuous density and pressure profiles, while for the CIB we approximate group and cluster galaxies as point-sources with a stochastic halo occupation distribution (HOD) model. The pipeline has been developed in such a way that other models, such as oriented anisotropic profiles that depend on the tidal tensor and more complicated dependence on the assembly history of halos and their environments, are straightforward to implement. A halo-based approach is not sufficient to include the large scale correlations induced by lensing and the kSZ effect, which have a non-negligible contribution from regions exterior to galaxy and cluster scale dark matter halos. We hereafter refer to all material outside of halos resolved by the simulation as the 'field' component. We use second-order Lagrangian perturbation theory to predict the electron momentum and total matter density of this field component, ensuring that large-scale kSZ and lensing correlations are included self-consistently.
We illustrate the separate treatment of simulation cells, or particles, interior and exterior to halos identified in this work with the mass-Peak Patch method in Figure 3.1. Material interior to halos is shown in red, and the field is shown in blue. The dynamics of both field and halo are treated by calculating displacements using Lagrangian perturbation theory. These displacements are used to move the field particles and halos to Eulerian space. Material inside halos is redistributed over either a diffuse radial profile (Section 3.1.1), as shown in red, or into distributions of subhalos (Section 3.1.2). The contrubution from each field particle is determined by the field redshift-volume kernel function described in Section 3.1.3.
N c e l l = Lagrangian Eulerian projection halo/field division Halo Field M a p = F h a l o + F fi e l d Figure 1. left -Illustration of the separate Lagrangian treatment of simulation cells (particles) interior and exterior to halos. Material interior to halos is shown in red, and the field is shown in blue, where the dynamics of both the field and halos are determined by LPT. Material inside halos is redistributed over either a diffuse radial profile (Section 3.1.1), as shown in red, or into distributions of subhalos (Section 3.1.2). The contribution from each field particle is determined by the field redshift-volume kernel function described in Section 3.1.3.
Full-sky component maps for the halo component are generated using the 'scatter' approach for the diffuse halo profiles, where each pixel is assumed to be smaller than the angular variation across the halo, while field and CIB maps are generated using a 'gather' approach, where the contribution from a given mass element or galaxy is treated as a point and therefore contributes to only a single pixel. In the case where either of these limits is not satisfied, pixels or lattice sites can be subdivided, depending on requirements, which is particularly easy to do with regular sky pixelizations such as Healpix [79,80] and the regular periodic lattices used to characterized the field displacements.
Generically, weak lensing convergence κ, kSZ temperature fluctuation ∆T /T , and tSZ Compton y-parameter y in a given directionn can each be expressed as a specific quantity, f , integrated along the line of sight, where W (z) is redshift-dependent kernel that encapsulates the physics associated with F. Since all the effects considered here are linear, the final map is a sum over individual components, for diffuse halo profiles, field, and point-sources, respectively. A schematic diagram of this separate treatment of halos and the field is shown in Figure 3.1.

Diffuse Halo Profiles
Once a halo light cone catalog has been created we use the linearity of equation 3.2 to optimize the generation of halo maps through the use of look-up tables. In the cases we present here, these tables specify the contribution from a halo with a given mass M, redshift z, and angular positionĥ, to the map at a given point in the skyn: , and f h (r/R h |M, z) encodes the complex astrophysics associated with, e.g., the matter or pressure distribution. All the diffuse halo profiles models adopted in this paper have the form of a generalized NFW profile (gNFW; [81]) with a mass and redshift dependence parameterized by . (3.4) For each halo i at angular positionĥ i in the catalog, its contribution to each pixel j is given by where the integration is over the solid angle subtended by the pixel, with δ Ω j ≡ Ω j dn j . We evaluate the values sampled at the pixel center to obtain an initial estimate of the value using and ensure the integrated contribution of the halo over all overlapping pixels is correct by an overall normalization factor: The total value for each pixel i in the halo map, F h i , is obtained by summing over all halos, The values of j for which δ F i j = 0 for a given i are computed using the condition We convert between peak patch halo masses and spherical overdensity masses fit in the hydrodynamical simulations by abundance matching along the light cone using the universal halo mass function of [82] with a mass 200 times the mean, M 200m . When necessary, we convert to halo masses with a different overdensity by assuming an NFW profile with a fixed concentration, c = 7, to be consistent with the mass profiles used in the lensing calculation (see §3.2.3).

Point Source Halo Occupation Distribution
We assume galaxies populating a given halo emit with the same spectral shape, separated into central and satellite galaxies. The distribution of subhalos of a given mass m, hosted in halos with mass M, is used to assign luminosities to satellites by assuming that each satellite's properties are determined by the subhalo mass and redshift. In particular we use the fit of subhalo mass function (SHMF) from N-body simulations in [83], . Note that this relationship only depends on the ratio of subhalo to host halo mass, and not on redshift or host halo mass, which is a good approximation given the considerable other uncertainties in modeling the properties of dusty infrared sources within halos. For each dark matter host halo in the light cone catalog, we draw from this distribution in two steps. First, we determine the mean number of subhalos contained in a halo of a given mass, Each halo is assumed to host one central infrared galaxy, following the infrared halo model discussed below in Section 3.2.2. The number of satellites in a halo with mass M is drawn from a Poisson distribution with mean of N sh (M), and its spatial distribution is chosen with a random angular position with respect to the halo center and with a radial distribution drawn from an NFW profile with a massredshift-concentration relationship given by [84]. Once the number of satellites has been sampled from the Poisson distribution, each one is assigned a subhalo mass, m, drawn randomly from the subhalo mass function consistent with equation (3.9) and a minimum mass of M min .

Anisotropies Generated Exterior to Resolved Halos
In order to include Thomson scattering and lensing by material external to halos, we use a simplified model for the spatial distribution of matter outside of halos that uses Lagrangian perturbation theory.
In particular, for each lattice site, i, in Lagrangian space that is outside of a halo at comoving position q i = |q i |n i , we determine the redshift corresponding to that point, z i = z(|q i |), assuming that a photon reaching the observer was emitted at that time from its Lagrangian redshift. As discussed in §2.4, the mass element belonging to that lattice site will have been displaced by peculiar motions, slightly changing the time at which the photon was emitted or scattered, but this is a negligible effect. Having determined the redshift corresponding to the mass element, we then determine its position as i and s (2) i are the first and second order LPT displacements as determined in §2.1. Once the perturbed position along the line of sight has been calculated, the total contribution from a given lattice site, δ F f i j , is determined by a redshift dependent kernel, such that where a 3 latt is the volume corresponding to each mass element or cell, v i = a i [Ḋ(z i )s (1) i +Ḋ (2) (z i )s (2) i ] is the peculiar velocity of the mass element, and Ω pix is the solid angle subtended by a pixel in the map. The factor W F (z, v) encapsulates the redshift of the source function and any line of sight peculiar velocity effects, and we describe the specific functional forms that we adopt for W F (z) in §3.2.4. We also include a generic Eulerian bias factor b F (z), suitable for modeling biased electron or dusty galaxy distributions, but have set b F = 1 in the lensing and kSZ maps presented here.
As lattice sites nearby the observer can subtend an angular extent larger than an individual pixel of the map, we calculate the angular size of each lattice site, and split the lattice site in three dimensions into the required number of sub-volumes (to a maximum of 5 3 for computational efficiency), such that the angular extent of a sub-volume becomes smaller than a pixel. The total contribution of the parent lattice site is then split evenly between the sub-volumes, and each sub-volume is independently added as a point to the map. This particle sub-division suppresses the shot-noise at small angular scales.

Astrophysical Models
Here we present the fiducial astrophysical models for determination of the frequency dependent CIB intensity, I cib ν (n), the thermal SZ Compton-y parameter, y(n), the kinetic SZ temperature fluctuation, ∆T ksz (n), and the CMB lensing convergence, κ cmb (n). These were determined from halo profile fits to cosmological simulations and empirical models from the literature. We note that, for example, the CIB has various halo models used in the literature [17,[85][86][87][88][89], but exploring the differences of these when applied to our full-sky simulations is beyond the scope of this work.

Thermal and Kinetic Sunyaev-Zel'dovich Effects
The thermal Sunyaev-Zel'dovich (tSZ) effect is determined by the Compton y-parameter, which is proportional to the optical depth-weighted ratio of electron thermal and rest-mass energies in a given direction: where P th is the thermal pressure of the gas, Y p 0.24 is the abundance of helium, and the medium is assumed to consist of fully ionized hydrogen and helium. This expression is valid up to small relativistic corrections when the electron temperature is much greater than the background radiation temperature, T e T γ . The spectral distortion induced has a characteristic effective temperature deviation given by [90] ∆T By using detailed hydrodynamical simulations including AGN feedback, [58] determined a parametric model for the pressure profile over a broader range of redshift and halo mass than is accessible by current observations. Their profile is in good agreement with the results from X-ray observations of nearby galaxy clusters [91], Bolocam observations of massive galaxy clusters [92], those from the Planck Collaboration [93], and stacked SZ profiles on locally brightest galaxies [94]. The dimensionless pressure profile is parameterized with the generalized NFW form given in equation , which follows the self-similar scaling for galaxy clusters [95]. The overall amplitude amplitude, F 0 , core-scale, x c , and large radius asymptotic power-law index, β , are fit for as functions of M, and z, while the other power-law parameters, α and γ, are fixed are fixed to 1 and -0.3, respectively. The pressure profile parameters are obtained from the values in Table 1 of [58], corresponding to AGN Feedback with ∆ = 200.
Doppler shifting of CMB photons scattered by free electrons along the line of sight results in a black body temperature fluctuation given by where τ is the Thomson scattering optical depth and n e is the comoving electron number density. The kSZ effect in our maps is separated into two components, a dense spherical halo and a clustered field, which includes electrons in the IGM, such that ∆T ksz (n) = ∆T ksz halo (n) + ∆T ksz field (n). To prevent double counting of collapsed matter in the field, we compensate the halo map by subtracting out the kSZ effect from a uniform sphere with the same position and mass as the halo, but with a radius corresponding to an overdensity of 3, somewhat smaller than than the Lagrangian radius of the halo, ensuring that mass and momentum is conserved on large scales.
We use halo gas density profiles fit to the same set of AGN simulations, with parameters given by [96]. The dimensionless gas density profile is also parameterized with the generalized NFW form given in equation (3.4), such that where ρ b is the mean comoving baryon density, and Y p = 0.24 is the helium mass fraction. Assuming all hydrogen and helium in gas is fully ionized, then F corresponds to the gas overdensity. Similar to the pressure profiles, the overall amplitude amplitude, F 0 , core-scale, x c , and large radius asymptotic power-law index, β , are fit for as functions of M, and z, while the other power-law parameters, α and γ, are fixed are fixed to 1 and -0.3, respectively. See §3.2.4 for the treatment of kSZ for the field component in material exterior to halos.

Cosmic Infrared Background
The CIB is produced by star-forming galaxies when stellar radiation is absorbed by dust grains and re-emitted in the infrared. Star formation is dependent on host halo mass, environment, and redshift, and is suppressed at low and high masses by supernovae, active galactic nuclei (AGN), and other kinds of feedback. Since the clustering of halos and the galaxies within them source the observed intensity fluctuations, the CIB provides important empirical constraints on the connection between star formation and redshift over wide range of halo masses and redshifts, in particular for the faintest and most difficult galaxies to detect and study individually. For this work we use the CIB halo model developed by [89] and used by [17] and [87]. In particular we adopt the parameters used by [87] to fit CIB power spectrum measurements with Herschel. We present the main details of the CIB halo model here, in addition to the elements that are unique to creating mock full-sky observations, but refer the reader to the literature for additional details and discussion [17,87,89]. In this model, the rest-frame SED of a given source depends on the frequency ν of observation, (sub)halo mass M, and the redshift z: where: • The spectral energy distribution, Θ[ν, T d ], is a greybody at low frequencies and a power law at high frequencies, Here B ν denotes the Planck function, and β = 1.6 is dependent on the physical nature of the dust. The relative normalization of the two components (since the absolute normalization is already accounted for in L 0 ) are obtained by the requirement that the logarithmic slope is continuous at ν 0 , dlnΘ(ν, z)/dlnν = −γ. Finally, the redshift dependencies of the effective dust temperature is given by T d ≡ T 0 (1+z) α , where T 0 = 20.7 and α = 0.2. For the simulations used in this work, z < 4.6, and the frequencies considered, ν < 1000, the spectral energy distribution remains a greybody.
• The redshift dependent global normalization of the L − M relation is of the form where δ CIB = 2.4. M e f f here describes the peak of the specific IR emissivity, and σ L/M describes the range of halo masses which produce the luminosity. We use the values of log(M eff /M ) = 12.3 and σ 2 L/M = 0.3 from model 1 of [87]. These act to produce a peak in the L − M relation describing the maximum average infrared emissivity per unit mass, due to the suppression at both the high and low mass end.
• L 0 is a free normalization parameter, which we chose to reproduce the Planck 2013 [17] power spectrum results at 545 GHz for = 500 in §4.4.3.
Having determined a population of subhalos for each host halo in the catalog, we populate each subhalo of mass m with a satellite galaxy with luminosity L ν (m, z), where z is the redshift of the parent halo. We also include the emission from a central galaxy for each halo of mass M at redshift z, and assign it a luminosity L ν (M, z). This treatment considers the L − M relation of satellites and centrals to be equivalent and ensures that every halo contains at least one infrared source.

CMB Lensing
As light propagates through the large scale structure of the universe it is deflected by the gravitational attraction of mass along the line of sight. In generating our maps we use the Born approximation, in which the displacement of light rays from their original position on the sky is small. We note that for upcoming surveys, post-Born terms will likely be non-negligible [97,98]; in principle, a detailed ray-tracing approach to CMB lensing that does not implement the Born approximation could be applied to our simulations, but we leave implementation of this to future work. Under the Born approximation, the effect of lensing on the primary CMB is completely specified by the convergence, κ, and we need only integrate the matter overdensity along the line of sight, weighted by the lensing kernel W κ (χ), where W κ (χ) = 3Ω M H 2 0 (1 + z)χ(1 − χ/χ * )/2 and χ * is the distance to the source. In this case, we choose the CMB as our source, located at a distance χ * 13.8 Gpc. The projection of matter density to generate convergence maps is analogous to the kSZ, except the latter only uses electron density, and is weighted by the velocity. Given the considerable astrophysical uncertainty in the shape of the matter density profile in the presence of feedback and star formation in the group and cluster mass halos resolved by CMB observations, the matter in all halos is assumed to follow an NFW profile with c ≡ r 200 /r s = 7, independent of mass and redshift, where r 200 is the radius enclosing 200 times the mean matter density. To account approximately for nonlinear infall in the outskirts of the halo, we extrapolate the profile using ρ(r) = ρ NFW (r 200 )(r/r 200 ) −2 for r 200 < r < 2r 200 . Additionally, halos whose virial radii subtend a solid angle less than twice that of a pixel in the map are considered unresolved and are included in the field component, along with all the other matter exterior to resolved halos. The distribution of field matter is determined as described in the next section. To prevent double counting, we follow a similar procedure that for the electron density in the kSZ calculation described in 3.2.1. This ensures that overall mass is conserved and the resulting convergence maps naturally obey κ = 0.
Material beyond the maximum redshift of the simulation (z = 4.5) is not resolved by either the halo or the field component, and is instead included by generating an uncorrelated gaussian random map with a power spectrum determined using the Limber approximation, and adding the z > 4.5 gaussian random map to the z < 4.5 halo+field map. We use the CAMB 'takahashi' Halofit model [99,100] with the parameter 'lens_potential_accuracy' set to 2 in order to calculate P m (k = /χ, z). We perform lensing of the CMB itself using the pixell package 4 . We first generate a Gaussian random CMB map of both temperature and polarization. We reproject the Websky convergence map in the Healpix pixelization to a lensing potential map in the cylindrical projection and pixellization used in the pixell package. We then perform the lensing, at 1' resolution. The statistics of the resulting lensed CMB are described in Sec. 4.4.4.

Scattering, and Lensing from Outside of Resolved Halos
As mentioned in §3.1.3, observables can originate from regions exterior to resolved halos. In particular, angular fluctuations on the largest scales in kinetic SZ and lensing are mostly determined by material outside of the group and cluster mass halos resolved in the simulations presented in this paper. Here we detail the specific expressions for going from the field matter distribution described in §3.1.3 to an observable map. In all cases, the effects are encapsulated by the field redshift-volume kernel function W f (z) defined by equation (3.11). These are: where f e accounts for the mean fraction of gas that is ionized with mean molecular weight of µ e . We account for the measurements of the cosmic mean electron fraction from hydrodynamical simulations in [101] in an approximate way by setting f e = 0.9. We assume hydrogen is fully ionized and helium is once ionized at z > 3 and fully ionized at z < 3 to set the mean molecular weight of ionized gas. In this work, electron and matter fluctuations in the field are assumed to follow Lagrangian perturbation theory, so b τ = b ksz = b κ = 1. Conversely, the halos hosting the CIB flux in our model are mostly resolved within our simulation, so we do not include a CIB field component.

Large-Scale Structure Simulation
The synthetic extra-galactic mocks presented in this work are constructed from a (15.4 Gpc) 3 , 12,288 3 particle lightcone generated using the mass-Peak Patch method [68,69], which we summarized in Section 2. This was accomplished by creating periodic initial conditions for a (7.7 Gpc) 3 , 6,144 3 particle simulation, placing an observer at each of the eight corners the volume, performing the eight runs in serial, and stitching the results together to generate a single seamless halo catalog over the full-sky. This 'octant' method has the advantage of doubling the one dimensional resolution of the simulation, or equivalently raising the volume resolution by a factor of 8, while using the same memory footprint. The disadvantage is that it replicates structures on scales equivalent to the box size, but although a given region in the volume is replicated eight times in this configuration, it will in general be observed at a different time from a different direction, making each octant of the sky unique, while avoiding any discontinuities or artefacts. Any periodic effects are limited to very large angular scales, < 10, and there is no repetition of structure along any given line of sight. The cosmology used for the simulations is consistent with the results of Planck 2018 [70]: Ω m = 0.31, Ω Λ = 0.69, Ω b = 0.049, h = 0.68, σ 8 = 0.81, and n s = 0.965, and the linear power spectrum used to create the initial density field was generated with CAMB 5 . The run time of the simulation was only 3.84 hours on 1128 Intel "Skylak" 2.4 GHz cores of SciNet's Niagara cluster [102], for a total runtime of 4336 hours and a peak memory footprint of 7.67 TB (∼ 9 floats per resolution element). 5.9TB of disk space was required to store the initial conditions for additional post-processing, and 33GB was required for the halo catalogue. The computational efficiency of the mass-Peak Patch method allows for many realizations of a simulation of this size, but this currently remains as future work.

Halo Catalogue
We consider all halos within a distance 7.7Gpc of the observer, corresponding to maximum redshift of z = 4.6 and a total survey volume of 1, 900 (Gpc) 3 . Although the mass function of mass-Peak Patch halo catalogues is generally within ∼ 10% of N-body [68], to ensure mass-matching with an equivalent N-body run we performed abundance matching along the lightcone to M 200ρ m of [82]. Abundance matching was performed by calculating the mass-Peak Patch mass function, N(> M|z), in redshift bins of ∆z = 0.1. For each halo mass bin, at each redshift, we determine the fractional change in mass needed for the number of mass-Peak Patch halos in the equivalent volume to match that predicted by [82], and save this in a two-dimensional table. The mass of each halo along the lightcone is then multiplied by a fractional change determined by bi-linear interpolation from this table. We retain halos in the catalogue with a pre-abundance matched mass greater than 10 particles of the simulation. This results in a redshift dependent mass completeness of the final catalogue of roughly M min,Websky 1.2 × 10 12 M between 0 < z < 4 and increasing to 4 × 10 12 M by z = 4.6. An exact tabulation of the halo mass resolution as a function of redshift is provided alongside the halo catalogue.
The final catalogue contains roughly 9 × 10 8  Shown in Figure 2 is a thin equatorial wedge of the halo catalogue subtending a degree, representing slightly less than a percent of all the halos in this realization. The cosmic web of groups and clusters can be easily identified in the right two panels, with filaments, sheets and voids several tens of Mpc across, and smaller groups clustered around clusters. The peculiar velocity from the observer, illustrated by the halo's colour, shows the coherent large-scale velocity flows on scales much larger than the density, of up to many hundreds of comoving Mpc. The octant method used to generate the simulation assures no repetition of structure or velocity information along the line of sight, and large-volume simulations such as the one presented here are required to realize accurate large-scale velocity flows.

Extragalactic Maps
Shown in Figure 3 are the full-sky maps determined by projecting the halo catalogue and initial density field using the four models described in §3. Clockwise from top left we show the thermal Sunyaev-Zel'dovich, kinetic Sunyaev-Zel'dovich, cosmic infrared background, and the CMB lensing convergence signals, with a 10 • × 10 • zoom in at the native resolution of the publicly available maps (N side = 4096) in the center. The correlation of different map components is clearly visible. Perhaps most prominent is the effect of nearby cluster-mass halos in the Compton-y map, illustrating the strong bias towards massive objects at low redshift that contribute most of the signal due to their very high temperatures and late formation times. The kSZ map also reflects electron density inhomogeneities associated with massive nearby clusters, on small and intermediate scales, with the overall effect arising from the coherent velocity flows that dominate the structure on much larger scales. Having a large simulation volume approaching the cosmic variance limit is essential for simulating the SZ sky, most importantly to sample the large scale velocity modes that contribute significantly to the kSZ effect and to probe the tail of the mass distribution of the rarest galaxy clusters that is most sensitively probed by the tSZ effect. The CIB and CMB lensing convergence, on the other hand, reflect structure from higher redshifts and lower halo masses, where both the infrared emission and number of galaxies per sky-area peak, and consequently those maps are more diffuse in appearance and appear more correlated with each other than they do with the tSZ or kSZ maps. Bright unresolved galaxies contribute a notable flux to individual pixels of the CIB map (white), while the CMB lensing convergence is more smooth on small scales due to the diffuse halo profiles and the unity bias of the matter density field. The CIB (convergence) map was smoothed using a 2 (4) arcminute beam for visualization purposes.
Using the Websky convergence map we lens a randomly generated primary CMB to create a lensed CMB map, as described in Section 3.2.3. This completes the set of extragalactic foreground maps presented in this work. To compare the extragalactic maps on an equal footing and illustrate their effects when added as foregrounds to the primordial CMB a few frequency dependent unit conversions are required, which we describe in Appendix A.  Figure 4 shows the full effects of the extragalactic foregrounds on observations of the primordial CMB at a range of Planck frequency channels from 100 GHz to 353 GHz. The mean temperature at each individual frequency has been subtracted, so that ∆T represents the fluctuations about the mean sky temperature at that frequency, not about T cmb . We show the primordial CMB in the top left panel, while the remaining panels on the left show the observed CMB after adding the four extragalactic foregrounds simulated in the Websky suite: tSZ, kSZ, CIB, and lensing. The top right panel shows the effect of lensing, while the other panels on the right hand side show the SZ and CIB foregrounds with no CMB contribution, using a different colour scale. Lensing was kept distinct to more easily distinguish the other foreground effects at low frequencies. The frequency dependent conversion factors for Compton-y and the CIB, coupled with the redshift evolving SED inherent in the CIB model, result in visually distinct effects in the observed CMB maps as a function of frequency.
Below 217 GHz, lensing is the dominant large-scale foreground, introducing a characteristic hot/cold variation wherever positive/negative features in the lensing convergence map happen to align with gradients in the primordial CMB. The standard deviation of the lensing effect is ∼20 µ K, with a minimum and maximum value pixel value of roughly ± 200. Focusing beyond lensing in the remaining panels on the right hand side, we see that the CIB is the dominant foreground at frequencies above 100GHz, as expected. By 353 GHz the CIB has begun to add significant signal to the primary CMB, and by 545 and 857 GHz the primary CMB and all other foregrounds are nearly indistinguishable from the CIB contribution. For this reason we do not show frequency channels above 353 GHz. The 100GHz foreground map has a contribution from both the CIB and the the tSZ effect -which at this frequency shows up as a temperature decrement due to the upscattering of photos to higher frequencies when passing through clusters. At 217 GHz the tSZ has nearly no net effect, as the amount of photons upscattering out of the infinitesimally narrow band used here is equal to the number upscattering into the band from smaller frequencies, and the CIB again dominates. Below 217 GHz the tSZ is most apparent, and at 100 GHz the large scale flow of the kSZ is also slightly visible, alongside a subdominant kSZ contribution from individual large halos with relatively large peculiar velocities.

Power Spectra
In this section we show the auto power spectra of thermal SZ, kinetic SZ, CMB lensing convergence, and the CIB at selected frequencies. We also show the cross-power spectra of the CIB with both thermal SZ and lensing, and the lensed TT, EE, and BB power spectra for a random primary CMB realization with the same cosmology. Where appropriate, we compare the power spectra to those from existing observations, simulations, and analytical predictions.  Figure 5. The power spectrum of the Websky Compton-y map (blue line). We also include the Planck measurement (black) [13], an external Planck re-analysis taking into account the trispectrum in the covariance matrix and an extended marginalization (gray) [103], the ACT measurement (green) [104], and the SPT measurement (orange) [105].

Thermal Sunyaev-Zel'dovich
Following the thermal Sunyaev-Zel'dovich halo prescription described in §3.2.1 we projected the pressure profiles determined from the hydrodynamical simulations of [58] onto the mass-Peak Patch halo catalogue to create the Websky Compton-y map. All ∼9 × 10 8 halos were used, and no contribution from the field was considered. Although the minimum halo mass used to determine the pressure profiles was ∼1×10 13 M , the uncertainty introduced by using these profiles for less massive Websky halos is neglegible compared to the total signal, due to the approximately M 5/3 scaling of pressure. Figure 5 shows the resulting power spectrum (blue line). For reference we include the Planck results (black) [13], and the external Planck re-analysis (gray) which took into account the trispectrum in the covariance matrix and an extended marginalization [103]. We also include the ACT measurement (green) of D ν=148GHz =3000 = 3.4 ± 1.4µK 2 [104], and the SPT measurement (orange) of D ν=143GHz =3000 = 4.08 +0. 58 −0.67 µK 2 [105]. We find good agreement with the three datasets at all values of . The large-scale tSZ power spectrum below ∼ 300 is highly dependent on the individual cosmological realization due to cosmic variance. In a similar study where we created ∼ 200 full-sky realizations of a smaller volume to z = 1.25 we found variance contours at < 100 similar in size to those shown for [103], but the Websky realization aligns well with Planck at these large scales. The small scale power spectrum ( > 10 3 ), on the other hand, is sensitive to the details of the halo pressure profiles used, as it directly encodes the shape of the gNFW profile, integrated over the redshift-evolving mass function and halo viewing angle. We find a slight increase in the power at = 3000 when compared to both the SPT and ACT measurements. Whether the disagreement is due to systematic effects in the SPT and ACT measurements, or indicates a physically relevant small-scale suppression in the tSZ power is not clear, due to the limited frequency coverage and sensitivity of SPT and ACT.
A free parameter when pasting on diffuse profiles is the integration radius at which to cutoff the halo profile. Due to the small radius power-law index of γ = −0.35, the pressure profile in the inner regions of halos falls more slowly with radius compared to the density profile, while beyond a radius of ∼ R ∆ the signal quickly falls due to large radius power-law index of β ∼ −4.35 (with a slight redshift dependence). Therefore, only integrating profiles to a small radius results in not accounting for signal from the outskirts of halos, and too large of radius will cause excessive doublecounting of signal from neighbouring halos. For example, when using R max = [2, 3, 4]R ∆ , the total Compton-y contribution of a 1 × 10 14 M halo at redshift zero is a factor of ∼ [1.9, 2.4, 2.7] larger than if stopping the integration at R max = R ∆ . In the power spectrum this effect mainly shows up as a slight normalization change for the range of multipoles shown, except for at the very smallest scales where the shape and extent of the profile is reflected. Since we do not attempt to model the thermal pressure distribution of material not associated with the resolved halos from the simulation, we chose an integration radius of 4R ∆ for this study. If one was to use the resulting Compton-y map for a detailed investigation of e.g. filamentary structures, the effects of this choice resulting in potential double counting of signal in the outskirts of halos would first need to be studied in further detail beyond the power spectrum analysis performed here.

Kinetic Sunyaev-Zel'dovich
The kinetic Sunyaev-Zel'dovich map contains both a spherical halo profile, and a clustered field component containing the electrons in the intergalactic medium and unresolved halos, described in §3.2.1 and §3.2.4, respectively. All halos for which r 200c subtends greater than half an arcminutue and M 200c is greater than 10 13 M were used in the creation of the halo map, and material in halos smaller than this was considered to be part of the field. This distinction was implemented due to the minimum halo mass used to determine the free electron profiles of M 200c ∼ ×10 13 M . The power spectra of the Websky halo (orange), field (blue), and total (gray) kSZ map are shown in Figure 6.  [101]. We depict the Battaglia et al. result with a different colour as it was measured from the same set of hydrodynamical simulations used to fit the free electron profiles applied in this work. The data points for all non-Websky simulations are adapted from Table  1 of [101], where we first used the radiative cooling and star formation (CSF) scaling relation given in Table 3 of [106] to scale their results to our cosmology, and assumed that the ratio of powers between models at = 3000 holds under the cosmological scaling performed. We find that the kSZ power spectrum is dominated at all scales by the contribution from the field component, with the halo component contributing roughly 5/6ths of the total signal at = 5000. The resulting total kSZ power spectrum is in rough agreement with the available simulation based predictions.
Given that the field component is based on second order Lagrangian perturbation theory, we expect the resulting power spectrum to be approximately consistent with the leading order contributions to the kSZ power spectrum from linear density and velocity fluctuations. The kSZ contribution in a given directionn can be written as  Table 1 of [101], where we first used the radiative cooling and star formation (CSF) scaling relation given in Table 3 of [106] to scale their results to our cosmology, and assumed that the ratio of powers between models at = 3000 holds under the cosmological scaling performed. We also include analytical calculations of the Doppler contribution (dashed black) resulting from the sharp cutoff of structure at z=4.5 in the Websky simulations, and the contribution from Ostriker-Vishniac component (solid black) [109].
δ e is the electron density contrast, δ is the gas density contrast, δ x is the ionization fraction contrast, v is the velocity, x(z) is the mean volume averaged ionized fraction, and the remaining terms are described in §3.2.1. The first three contributions are commonly referred to as the 'Doppler' (∝ v), the 'Ostriker-Vishniac' (∝ vδ ), and the 'patchy' (∝ vδ x ), respectively [109]. Although reionization is not included in these z<4.5 simulations, we find a large scale power of a few µK 2 resulting from the sharp cutoff in structure at the maximum redshift of the simulations, z max . This sharp cutoff suppresses cancellation of the line of sight velocities from comoving distances beyond the bounds of the simulation, with a power given by where u 0 (z) ≡ σ T n e,0 (1 + z)H f and f ≡ d ln D/d ln a. When integrating to a maximum redshift of z max = 4.5, this is the dominant source of power at large angular scales. The Ostriker-Vishniac power spectrum dominates on small scales and is given by where µ ≡k ·k . We include the Doppler and Ostriker-Vishniac components in Figure 6. Previous Nbody simulations of the kSZ sky at z<1 [46] contained similar descriptions of the field (model I), halo (model II), and total (model III) as shown here. While the different redshift coverage, model for the kSZ signal of clusters, radial extent of the cluster profile integration, and cosmological parameters, make a direct comparison impossible, we find relatively good agreement.

Cosmic Infrared Background
The Cosmic Infrared Background is simulated using the halo occupation distribution model described in §3.2.2. We populated all ∼9 × 10 8 halos in the mass-Peak Patch catalogue with a number of central and satellite galaxies, resulting in a total of ∼ 3 × 10 9 galaxies projected into the final map. Due to the mass resolution of the Websky halo catalogue of ∼ 1 − 2 × 10 12 h −1 M , the CIB contribution from the field was not considered for this study. Maps were produced at the Planck HFI frequencies of 100, 143, 217, 353, 545, and 853 GHz, and a number of additional frequencies corresponding to the bands of current and future ground-based CMB experiments around this range. The Planck CIB model used requires the fitting of an overall normalization factor L 0 , which we determined by requiring the power spectrum of the 545 GHz Websky map to match the 545 GHz CIB auto-power spectrum at = 500, as measured by the Planck Collaboration [17]. Figure 7 shows the power spectrum of the Websky CIB maps at the Planck frequencies with available data: 143, 217, 353, 545, and 853 GHz. We also include the 2013 results of the Planck collaboration (black circles) with the assumed radio-source shot-noise levels of 6.05 ± 1.47, 3.12 ± 0.79, 3.28 ± 0.82, 2.86 ± 0.7, and 4.28 ± 0.90 subtracted out [17], the Herschel Multi-tiered Extragalactic Survey (HerMES) measurements of Viero et al. (green triangles) [87], and the recent dust-cleaned Planck maps of Lenz et al. (gray pentagons) [110]. The HerMES results were obtained using the SPIRE instrument at 250, 350, and 500 µm (1200, 857, and 600 GHz). Here we include the 857 GHz results as downward-facing triangles, and the 600 GHz results (note not 545 GHz) as upward-facing triangles, both using the 'only extended sources masked' measurements provided. The Lenz Table 1. Frequency decoherence of the CIB measured by averaging C νν /(C νν C ν ν ) 1/2 over the range 150 < < 1000. Error bars correspond to the standard deviation in this range. We include the Planck measurements of [17] and the Lenz et al. measurements of [110].
Planck analysis used additional data of neutral atomic hydrogen from the recently-released HI4PI Survey to create and remove template maps of Galactic dust, allowing for a measurements to extend to larger scales. We show the baseline results of Lenz et al., which used a low H I column density mask corresponding to a threshold of N HI = 2.5 × 10 20 cm −2 . We do not show their results at > 700 as they are effectively the same as those from the Planck at small angular scales.
We find good agreement on the shape and amplitude of the CIB power spectrum at all frequencies, with a slight excess of power on the Poisson tail out to = 10, 000 at 857 GHz. The Planck and Hermes data shown here do not have equal flux density cuts applied. HerMes 'extended sources masked' corresponds to a frequency-independent flux density cut of 400 mJy, while Planck uses frequency dependent flux density cuts of [350,225,315,350, 710] mJy for [143, 217, 353, 545, and 857] GHz, respectively. For the purpose of the CIB power spectrum plots shown in this paper, we applied a simplistic 400 mJy flux density cut to the nside=4096 Websky maps, by setting the value of any pixel with a flux density greater than 400 mJy to the mean of the remaining pixels below the flux density cut. This resulted in nulling [3,109,1229,12772, 96011] pixels out of 200 million for the [143, 217, 353, 545, and 857] GHz maps. Additionally, the Websky data shown in this work uses the ideal case of an infinitely narrow passband centered at each frequency, while Planck and HerMes inherently include different spectral response functions integrated over the width of the passband and some assumption about the shape of the source spectrum. Although not included here, this can have a ∼30% effect on the power spectrum between the two observations [87], and one can expect a similar effect on the Websky results when treated equivalently.
As well as the powersectrum measurements, we determine the frequency decoherence of the Websky CIB maps through the cross correlation C νν /(C νν C ν ν ) 1/2 . The CIB is expected to decohere between channels due to the redshift-dependent peak in the CIB distribution. We report this quantity averaged over the range 150 < < 1000 in Table 1, where we restrict to this range to focus on the clustered CIB contribution and negate the contribution from shot noise as done in the Planck 2013 analysis. We neglect the decoherence with the 143 GHz map due to the high shot-noise contribution dominated by the radio sources [17] which were not included in our study. We find similar results to the Planck measurements of [17] and the Lenz et al. measurements of [110], with significant decoherence between 857 and 217 GHz, and a small but non-zero decoherence for neighbouring Figure 8. The lensing convergence power spectrum of the total Websky map (gray) and the Halofit model with cosmic variance errorbars (red). We also show the contribution from the three components necessary to construct the Websky map: halos (orange), field (blue), and the z > 4.5 contribution from material beyond the redshift resolved in the dark matter simulations used (green). frequencies.

Weak Gravitational Lensing
The CMB lensing convergence map κ and its lensing of the primary CMB are described in 3.2.3. The Websky CMB convergence map is constructed using the field particles and an NFW profile for all halos in the mass-Peak Patch halo catalogue that subtend more than a pixel, since lower mass halos are not resolved. We have added an additional component to the convergence, as an uncorrelated Gaussian random field, to account for the power originating from redshifts larger than the maximum included in the simulation, z max =4.5.
Given the matter power spectrum as a function of redshift, P m (k, z) the angular power spectrum of the lensing convergence, κ, is well-approximated by the Limber approximation of Equation 3.21. Shown in Figure 8 is the convergence found from Equation 3.21, using the 'Halofit' model for the matter power spectrum [99,100] integrated from z = 0 to z ≈ 1100 (solid red), and integrated from z = 4.5 to z ≈ 1100 (dashed red). Also shown are the power spectra from the Websky halo (orange), field (blue), and gaussian z > 4.5 (green) contributions, as well as the sum of the three (gray). The Websky and Halofit convergence power spectrum are consistent to within a few percent at < 1000, with a residual suppression reaching ∼ 20 percent at the smallest scales. This level of deviation from halo fit would be expected from 10 percent variations in halo density profiles that could be due, for example, to baryonic feedback effects. The large scale power is nearly entirely coming from material outside of resolved halos (field), while at > 3000 the halo contribution to the power dominates. As expected, the field power spectrum follows closely that of the linear matter power spectrum on large scales.
Using the CMB convergence map, we lensed a randomly generated primary CMB following the procedure described in section 3.2.3, and show the power spectrum results of the simulations (top) and comparison to theory (bottom) in Figure 9. The simulation results for the temperature (TT) and . CMB lensing power spectra. Shown in the top panel are TT, EE and BB of unlensed (solid) and lensed (dashed) primary CMB temperature and polarization anisotropies, based on the maps generated with the procedure described in section 3.2.3. The unlensed BB power is zero since the unlensed primary maps were generated using a transfer function with tensor-to-scalar ratio r = 0. Clearly evident is the well-known 'peaksmearing' effect lensing has on the acoustic peaks, and the generation of lensing B-modes. The lower-left panel shows the ratio of lensed to unlensed power for TT and EE for our maps (dots), along with the prediction from CAMB (dashed and solid lines for EE and TT, respectively), with the same Halofit model for the matter power spectrum that was used in Figure 8. The lower-right panel shows the same BB power from the lensed CMB map as shown in the top panel (dots), as well as the prediction from CAMB and Halofit. We find excellent agreement at all < 5000.
E-mode (EE) polarization show the well-known peak-smearing effect of lensing, where the acoustic peaks of the primary CMB (solid) become washed out into the lensed CMB (dashed), due to the local magnification and de-magnification effects of lensing which acts to transfer power between scales. The B-mode (BB) power of the unlensed CMB is identically zero due to the scalar-to-tensor ratio used as input, r = 0, but the gravitational-lensing generated B-mode component is non-zero and matches well with theory, as seen in the bottom right panel. The bottom left panel shows the ratio of the lensed to unlensed power for both TT and EE for the simulations (coloured dots) and the theoretical calculation from CAMB (solid and dashed line for TT and EE, respectively), using the same Halofit model as for the kappa analysis discussed directly above. We find excellent agreement for TT, EE, and BB, when compared to the CAMB results for < 5000. 100 GHz Figure 10. Cross correlation of the CIB and the CMB lensing potential φ , including the Websky results (orange), the Planck 2013 measurements (black circles) [17], and the recent large-scale CIB measurements from dust-cleaned Planck maps of [110] cross correlated with the Planck lensing convergence (gray pentagons).

CIB-Lensing Cross-correlation
The contribution of matter fluctuations of a given size to the CMB gravitational lensing potential φ , related to the lensing convergence κ through peaks about half way (z ≈ 2) between us and the last scattering surface. Similarly, the CIB redshift distribution peaks approximately around z 2 (with a dependence on observation frequency, being biased towards higher redshift at lower frequencies), and the dusty star-forming galaxies comprising it tend to reside in halos of mass 10 10 < M/M < 10 13 . Therefore, a significant correlation is expected between the two maps [111,112]. Observationally, detecting this correlation can be used to constrain models of the CIB, while here we use it to show that the Websky simulations properly encapsulate the correlation between maps necessary to provide accurate mocks for CMB analyses. This is not guaranteed from the power spectrum results of the Websky CIB and lensing convergence independently. The Figure 10 shows the cross correlation of the Websky CIB and gravitational potential maps (orange), the Planck 2013 results (black) of [17] which used the gravitational potential reconstructed in [113], and the recent large-scale CIB measurements from dust-cleaned Planck maps of [110] cross correlated with the Planck lensing convergence (gray pentagons). Due to the approximately −1 and −2 shape of the CIB and gravitational potential power spectra, respectively, the cross correlation C CIB×φ is typically plotted with a leading factor of 3 , as is shown here. We find good agreement at all frequency channels for the range of measured by Planck, which did not consider modes with < 100 or > 2000 in their analysis due to possible mean-field systematic effects and extragalactic foreground contamination, respectively. Correlated noise between channels can enter in a number of ways from, e.g., foregrounds in multiple channels, meaning that the measurements and errorbars at each frequency channel are not completely independent, as evidenced by first two Planck data points at low frequencies. We also show the recent dust-cleaned Planck maps of Lenz et al. [110]. Again we show the baseline results, which used a low H I column density mask corresponding to a threshold of N HI = 2.5 × 10 20 cm −2 . We find good agreement for both the amplitude and shape of the Websky  Figure 11. Cross correlation of the tSZ and 'sky frequency maps' T from 143 to 857 GHz in units of Comptony, where a frequency map is constructed by adding the CIB and Compton-y. We show the Websky results (orange) and Planck 2015 CIB-cleaned results [86] (black circles). The Planck uncertainties are dominated by foreground residuals, and thus they are highly correlated from one channel to another. and Lens et al. cross correlation, extending far below the original minimum Planck cutoff at these frequencies of 400. The turn-over at low also matches very well, which was not possible to determine from the Planck analysis alone.

CIB-tSZ Cross-correlation
The tSZ effect, unlike the CIB and gravitational lensing, is mainly produced by local high mass clusters. Therefore, the tSZ and CIB should have a much smaller overlap in redshift and halo mass than the CIB-lensing presented in the previous section. This small overlap makes the correlation hard to detect in practice, but determining the strength of the correlation remains important to determine the CIB contamination of tSZ estimates, to estimate the power resulting from the kSZ, and to constrain halo model approaches. In Figure 11 we show the CIB-tSZ cross correlation of the Websky maps (coloured) and the Planck 2015 results of [86] constructed by performing a cross spectrum with the Planck frequency maps and a reconstructed MILCA y-map derived from component separation, and subtracting out the CIB leakage estimated from simulations (see [86] for more details). The Planck data at 217 GHz was not available. The intensity units for the cross correlation, C yT ν , are "Compton y-units" [8], defined for each frequency bandpass i such that with conversion factors taken from Table 6 of [114]. See [86] for more details on how the crosscorrelation was calculated. We find broad agreement at all frequencies and multipoles for which there is significant signal-to-noise in the Planck measurements. At low frequencies the cross correlation is nearly completely determined by the Compton-y component due to the relatively low amplitude of the CIB, and we see the same turn-over at high as the tSZ power spectrum presented in Figure 5. By 353 GHz the CIB component begins to contribute roughly equal power as the tSZ, and by 545 and 857 GHz it dominates. See [86] for a discussion of limitations due to systematic effects and the statistical significance of the cross-correlation as measured by the Planck Collaboration.

Discussion
In this paper we have described a new pipeline for generating simulated maps of the microwave sky and catalogs of galaxy groups and clusters. Because it explicitly models the fate of individual connected 'patches' in Lagrangian space as a function of local measurements involving the strain of the linear displacement field, it is particularly well suited to connect observables efficiently with the initial conditions, without resorting to more costly N-body simulations. The pipeline consists of two main steps: first generate an accurate realization of large-scale structure on the past light cone and then use this realization to determine the observed spectral intensity and lensing convergence for each pixel in the simulated map. While the pipeline is applicable to a broad range of redshifts and scales, we have focused on a particular large-volume light-cone realization constructed by placing the observer at the eight corners of a periodic volume (∼ 5.25 Gpc/h on a side with N = 6, 144 3 resolution elements), resulting in a single realization with the observer at the center of a cube ∼ 10.5 Gpc/h on a side with N = 12, 288 3 resolution elements, from which the inscribed sphere was used to generate full-sky maps and halo catalogs out to a redshift z = 4.6. This simulation resolved halos above a mass of approximately 10 12 M , sufficient to include the complete tSZ signal without considering a contribution from the field component, although some signal from the low-mass end of the CIB, which is generally considered to have a minimum contributing halo mass of 10 10 -10 11 M , is missed. The free normalization parameter in the Planck CIB model used allowed for the Websky CIB maps to nevertheless reproduce the proper clustering statistics, but the omission of these low mass halos will change the average CIB bias, and is a likely reason that the Websky CIB power spectrum does not match the Planck results as well as a model that includes halos with mass less than ∼ 10 12 M /h. This signal could be included through the additional use of a field prescription, but was not included in this work. We also note that the inclusion of radio sources for the same halo catalog presented here is the subject of work in preparation.
Spherically symmetric halo profiles were used in all the maps presented here. Since in reality individual groups and clusters are not spherical, this is a limitation that should be kept in mind when using these maps. While the halo catalogs upon which the maps are based do not contain direct information about the final state (Eulerian) orientation, the measurements performed in the mass-Peak Patch method already contain a wealth of information of Lagrangian measurements around each halo on a number of scales. This information can be used to estimate of the final shape, orientation, and internal properties, and is a natural extension of our approach. Although we assumed the interiors of halos to be spherical, the halo clustering and associated matter and thermal energy distribution anisotropies in the two-halo regime are accurately represented in these simulations.
Additional tests of the mass-Peak Patch simulation method we used to construct the halo catalogues have been shown to accurately reproduce higher order halo statistics such as the halo bispectrum [77]. Higher order spectra of the Websky CIB maps have also been presented in [115], while correlations between the halo distribution and both CIB and Compton-y maps were shown in [116]. The map power-spectrum results and cross-spectrum results were compared to available data from Planck [13,17,86], Planck external re-analyses [103,110], ACT [104], SPT [105], the Herschel Multi-tiered Extragalactic Survey (HerMES) [87], and a number of hydrodynamical simulation results [43,52,101,[106][107][108] with broad consistency found for the following validation statistics: power spectra of tSZ, kSZ, CIB, and the lensing convergence; the CIB -CMB lensing cross-correlation; and the CIB -tSZ cross-correlation. The halo catalogue and maps presented here are publicly available online 7 and are included as part of the so_pysm_models framework 8 Table 2. Frequency dependent conversion factors for CIB maps in MJy sr −1 , and Compton-y maps, to a deviation ∆T [µK] from the primordial CMB blackbody of T cmb = 2.7255 K. These conversion factors assume an infinitely narrow passband centered at the given frequency. Multiply a CIB or Compton-y map by the given value to convert it.