CMB-S4: Forecasting Constraints on Primordial Gravitational Waves

Determining the origin of structure in the universe is one of the most important open problems in cosmology. CMB anisotropies sourced by early-universe density perturbations are currently the most powerful observational probe of the earliest mechanisms of structure formation. It is possible that the same processes that produced the density perturbations also sourced tensor perturbations, or primordial gravitational waves (PGWs). If this is the case, detecting a PGW signal would yield insight into physics far earlier than the epoch of recombination, and allow us to build an unprecedented understanding of the earliest moments of time.

PGWs leave imprints on the polarization of the CMB (Hu & White 1997; Kamionkowski et al. 1997; Seljak & Zaldarriaga 1997). In particular, the sensitivity of CMB measurements to gravitational waves arises from the generation of polarization at the surface of last scattering: to first order, scalar perturbations produce only even-parity E-mode polarization, while tensor perturbations produce odd-parity B-mode polarization as well. Thus, a measurement of primordial B-mode polarization in the CMB, parameterized by the tensor-to-scalar ratio r, is a direct measurement of the amplitude of tensor perturbations. Detecting r has profound implications for high-energy physics and the quantum nature of gravity (Krauss & Wilczek 2014), and the potential to shed light on the mechanism that produced these primordial perturbations.

Cosmic inflation is our current leading paradigm for what occurred in the very early universe. It was first put forward to explain the lack of observed magnetic monopoles and to solve the flatness and horizon problems (Kazanas 1980; Starobinsky 1980; Sato 1981; Guth 1981; Linde 1982, 1983; Albrecht & Steinhardt 1982) and has since been an active field of research. The theory describes a period of exponential expansion in which quantum fluctuations are magnified to cosmic size and become the seeds for all structure in the universe (Mukhanov & Chibisov 1981, 1982; Guth & Pi 1982; Hawking 1982; Starobinsky 1982; Mukhanov 1985; Bardeen et al. 1986). In addition to the production of PGWs (for a recent review, see Kamionkowski & Kovetz 2016), inflation makes several predictions, most of which—superhorizon fluctuations, Gaussian perturbations, adiabatic fluctuations, spatial flatness, and a nearly scale-invariant scalar spectral tilt—have been confirmed, most recently by the Planck collaboration (Planck Collaboration et al. 2020b).

There are currently a number of ground-based experiments measuring the CMB polarization to high precision on a range of scales, and attempting to constrain the tensor-to-scalar ratio, including ACT (Aiola et al. 2020), BICEP/Keck (BICEP2/Keck Array Collaborations X 2018), CLASS (Harrington et al. 2016), POLARBEAR/Simons Array (Suzuki et al. 2016; Hasegawa et al. 2018), and SPT (Bender et al. 2018; Sayre et al. 2020), with Simons Observatory to follow soon (Ade et al. 2019). Additionally, there are current and future balloon and satellite missions such as SPIDER (Gualtieri et al. 2018) and LiteBIRD (Hazumi et al. 2019), which we expect to complement ground-based measurements. The current best constraints are r(k_s = 0.05 Mpc⁻¹, n_t = 0) < 0.06 at 95% confidence (BICEP2/Keck Array Collaborations X 2018), where k_s is the scalar pivot scale and n_t is the spectral index of the tensor modes. Ongoing efforts in both the Atacama Desert and at the South Pole, between now and the start of CMB-S4, will steadily improve these constraints while continuing to prove the methodologies on which CMB-S4 will rely.

CMB-S4, anticipated to start observations in 2027, is intended to be the definitive ground-based CMB polarization experiment. It is designed to cross critical thresholds in constraining the B-mode polarization signature of primordial gravitational waves and in sensitivity to new light relics, while also improving our understanding of the nature of dark energy and General Relativity on large scales (Abazajian et al. 2016). To achieve these goals requires a significant increase in sensitivity, from 2 to 4 × 10⁴ detectors in Stage 3 experiments to roughly 5 × 10⁵ detectors, and an unparalleled control over other sources of signal (e.g., Galactic foregrounds, gravitational lensing, etc.) and of systematics. Therefore, CMB-S4 will require telescopes at multiple frequencies, each with a maximally outfitted focal plane of pixels utilizing superconducting, photon-noise-limited detectors, and likely novel analysis techniques. To understand the optimal design for achieving the desired science goals, in particular focusing on primordial gravitational waves, we present the development of the CMB-S4 r forecasting framework and its application toward determining the CMB-S4 baseline r survey.

The sensitivity achieved by a CMB experiment, which observes for a given number of years with a given number of detectors, is subject to a number of efficiency factors. These include the fraction of detectors that are actually functional, the achieved sensitivity per detector versus model prediction, the fraction of days per year spent observing, the fraction of observing time spent "on field," and the fraction of data passing weather and other cuts. As an illustrative example, if each of the above efficiency factors were 0.65 then the product is 0.65⁵ ∼ 0.1. Forecasts for the sensitivity to r, where excess low-frequency noise or systematic contamination can lead to additional (potentially large) sensitivity loss, are particularly challenging. Most previous forecasts of the sensitivity of CMB experiments to r have been ab initio—requiring assumptions as to the value of all efficiency factors. Possibly due to the human tendency toward optimism, many of these forecasts have not in practice been achieved once the data have been taken and analyzed: ABS (Essinger-Hileman 2011; Kusaka et al. 2018); BICEP/Keck (Hivon & Kamionkowski 2002; Yoon et al. 2006; Nguyen et al. 2008; BICEP2/Keck Array Collaborations VI 2016); EBEX (Reichborn-Kjennerud et al. 2010; Abitbol et al. 2018); Planck (Planck Collaboration 2006; Planck Collaboration et al. 2020b); POLARBEAR (Lee et al. 2008; Kermish et al. 2012; Polarbear Collaboration et al. 2020); QUIET (Lawrence et al. 2004; Araujo et al. 2012); SPIDER (Montroy et al. 2006; Fraisse et al. 2013; Gualtieri et al. 2018); SPTpol (McMahon et al. 2009; Austermann et al. 2012; Sayre et al. 2020). Therefore, for CMB-S4, we take an alternate approach, scaling from the overall performance achieved in the best available existing experimental results. We scale directly from published B-mode noise spectra and bandpower covariance matrices derived from multiyear maps that have passed systematics null tests. This bypasses the need for an unbiased, individual accounting of the various efficiency factors, and naturally incorporates all effects that impact real-world experiments. This also differentiates our current approach from forecasts other groups have made for CMB-S4 prior to the formal existence of our collaboration (Wu et al. 2014; Errard et al.2016; Barron et al. 2018).

1.1. The CMB-S4 r Forecasting Loop

The CMB-S4 baseline r survey, described below, has been continuously evolving by translating science requirements to measurement and instrument requirements, based on our understanding of the impact of astrophysical foregrounds, instrumental systematics, delensing nonidealities, and analysis methodology. To maintain forecasting realism as complexity increases, our general approach has been an iterative one. We rely on a closed forecasting loop, presented in Figure 1, to tie the semianalytic tools (which allow for fast optimizations) to map-based studies (which can include multiple layers of additional complexity). To establish our measurement requirements and the baseline experiment configurations that can achieve them, we perform multiple passes through this loop. In the figure, the section number accompanying each stage of the loop indicates the paper section in which that stage is described in detail.

Zoom In Zoom Out Reset image size

The main steps describing this process are as follows.

• 1.  
  Develop a semianalytic power-spectrum-level forecast, assuming noise performance that is scaled from analyses of real experiments.
• 2.  
  Use this forecasting tool to optimize the allocation of detector effort across observing frequencies, determining certain "checkpoints" in survey definition space.
• 3.  
  Use the checkpoint configurations to create standardized, version-numbered, map-based data challenges (DCs) for validation.
• 4.  
  Estimate science parameters from the DC maps with independent component-separation analysis methods.
• 5.  
  Check that independent analyses show recovery of science parameters from these challenge maps that match analytic forecasts, either in terms of variance or bias. If they do not, we revise the forecasts accordingly.
• 6.  
  Iterate steps 1–5, injecting increasing realism in the form of: (a) sky model complexity informed by the latest data and modeling efforts; (b) survey coverage based on proven observing strategies; (c) systematics whose form, parameterization, and likely amplitude are likewise guided by real-world experience; and (d) treatment of lensing.

1.2. CMB-S4 r Forecasting Workflow and Evolution

For the CMB-S4 Science Book (Abazajian et al. 2016), we developed a semianalytic forecasting framework specifically targeted toward optimizing sensitivity to the tensor-to-scalar ratio, r, in the presence of Galactic foregrounds and gravitational lensing of the CMB. Currently, only small-aperture telescopes (SATs) have reached the level of systematics control and noise performance necessary to pursue a ground-based, high-precision measurement of B-mode polarization down to low multipoles (ℓ ≃ 30), targeting the ℓ ≃ 80 peak from the polarization signature generated by PGWs at the epoch of recombination. Therefore, to forecast the performance of next-generation SATs, this machinery is based on scaling the bandpower covariance matrices (BPCMs) and noise power spectra (N_ℓ ) of current published SAT analyses, such as those from BICEP/Keck (BICEP2/Keck Array Collaborations VI 2016; BICEP2/Keck Array Collaborations X 2018). This automatically builds into the forecast all real-world inefficiencies including (but not limited to): imperfect detector yield; nonuniform detector performance; read-out noise; observing inefficiency; losses due to timestream filtering; beam smoothing; and nonuniform sky coverage.

At its core, this code is based on the BICEP/Keck parametric power-spectrum-based likelihood analysis. Such types of parametric analyses have also been extensively used by the Planck collaboration for ℓ ≥ 50 (Planck Collaboration XV 2014; Planck Collaboration XI 2016; Planck Collaboration et al. 2020a). We validate this approach using map-based simulations as part of the forecasting loop, and present detailed results in Section 6. Our confidence in the projections is grounded in the connection to achieved performance and published results.

Figure 2 presents a schematic representation of the semianalytic forecasting framework, identifying the user inputs, code modules, and outputs. The subsections that follow describe this framework in detail.

Zoom In Zoom Out Reset image size

2.1. Fisher Formalism

Given a likelihood function of the form

$$\begin{eqnarray}&&L(\theta ;d)\propto \displaystyle \frac{\exp \left[-\tfrac{1}{2}{\left(d-\mu (\theta )\right)}^{{\mathsf{T}}}{\rm{\Sigma }}{\left(\theta \right)}^{-1}(d-\mu (\theta ))\right]}{\sqrt{\det ({\rm{\Sigma }}(\theta ))}},\end{eqnarray} \tag{ 1 }$$

where d are the data bandpowers, θ are the theory parameters, μ(θ) and Σ(θ) are the bandpower expectation values and the bandpower covariance matrix given the parameters, and T denotes the transpose, we can calculate the expectation value of the log-likelihood curvature, evaluated at the position of the best-fit model:

$$\begin{eqnarray}&&{F}_{{ij}}=-\left\langle \displaystyle \frac{{\partial }^{2}\mathrm{log}L(\theta ;d)}{\partial {\theta }_{i}\partial {\theta }_{j}}\right\rangle .\end{eqnarray} \tag{ 2 }$$

This quantity, called the Fisher information matrix, measures how steeply the likelihood falls as we move away from the best-fit model, and F⁻¹ can be thought of as the best possible covariance matrix for the measurement errors on the parameters θ_i . It can be shown that $\sqrt{{({F}^{-1})}_{{ii}}}$ is the minimum obtainable standard deviation on the desired parameters (e.g., Cramér 1946; Kendall & Stuart 1979; Tegmark et al. 1997).

Inserting Equation (1) into Equation (2) yields

$$\begin{eqnarray}&&{F}_{{ij}}=\displaystyle \frac{\partial {\mu }^{{\mathsf{T}}}}{\partial {\theta }_{i}}{{\rm{\Sigma }}}^{-1}\displaystyle \frac{\partial \mu }{\partial {\theta }_{j}}+\displaystyle \frac{1}{2}\mathrm{Tr}\left({{\rm{\Sigma }}}^{-1}\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial {\theta }_{i}}{{\rm{\Sigma }}}^{-1}\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial {\theta }_{j}}\right).\end{eqnarray} \tag{ 3 }$$

We then calculate our parameter constraints as

$$\begin{eqnarray}&&{\sigma }_{i}=\sqrt{{\left({F}^{-1}\right)}_{{ii}}}.\end{eqnarray} \tag{ 4 }$$

In all the projections below, for each fiducial model considered, we choose to fix the covariance matrix with respect to the theory parameters, i.e., Σ(θ) = Σ, making the second term of Equation (3) identically zero. ⁹⁶ Equation (3) provides a clear picture of how the construction of the covariance matrix Σ impacts the final constraints, and how its misestimation could lead to constraints that are far too optimistic. It is with this in mind that we have decided to compute our BPCMs by directly scaling the values in achieved covariance matrices (see Section 2.5).

2.2. Forecasting Inputs

2.3. CMB-S4 Scalable Instrument Specification

To span the four atmospheric windows available to ground-based experiments (Figure 3) and have enough channels to mitigate against complex foregrounds, we assume eight channels at 30, 40, 85, 95, 145, 155, 220, and 270 GHz, which are placed on small-aperture telescopes. For these low-resolution instruments, we pick 0.52 m apertures, motivated by proven SAT Stage-3 experiments (Kang et al. 2018), yielding the beams presented in Table 1. In addition, we also include a 20 GHz channel on a large-aperture telescope (LAT). ⁹⁷

Zoom In Zoom Out Reset image size

Table 1. Optimized Instrument Configuration for the PGW Survey, as Presented in the CMB-S4 CDT Report

		Frequency (GHz)
Science Goal	Item	20	30	40	85	95	145	155	220	270	DL	Total
r	No. of detectors	130	260	470	17 k	21 k	18 k	21 k	34 k	54 k	84 k	250 k
	Angular resolution [FWHM]	11'	77'	58'	27'	24'	16'	15'	11'	85	10

Download table as:  ASCIITypeset image

We derive the split in each window by separating the overlapping bands as far as possible while still keeping the calculated per-detector NET within 10%–15% of the NET for a detector that spans the full window. The ideal per-detector NETs were calculated with NETlib.py ^{98 , 99} at the South Pole and Chile, using the 10 yr MERRA2 median atmospheric profiles (Gelaro et al. 2017). We use the average of NET calculated for detectors at the two sites, which are 214, 177, 224, 270, 238, 309, 331, 747, and 1281 μK_CMB $\sqrt{{\rm{s}}}$ for our nine channels at 20–270 GHz, respectively. ¹⁰⁰ These NETs are lower than similar BICEP/Keck channels (which are 287.6, 313.1, and 837.7 μK_CMB $\sqrt{{\rm{s}}}$ for 95, 150, and 220 GHz, predominantly because they are calculated for a 100 mK thermal bath, as opposed to 250 mK (which was also used in the CMB-S4 Science Book forecasts).

We also fold in information from two WMAP channels, 23 and 33 GHz, and seven Planck channels, 30, 44, 70, 100, 143, 217, and 353 GHz, though this extra information is only relevant in the early stages of CMB-S4 observation.

Since the forecasts use scaled BICEP/Keck bandpower statistics, they also use the same bandpower window functions and uniform binning: nine multipole bins with Δℓ = 35, spanning a multipole range of 21 ≤ ℓ ≤ 335.

2.3.1. Delensing Treatment

One of the main challenges for detecting primordial B modes is the lensing B-mode contribution from the weak lensing of E modes as the CMB photons travel to us. For instrument noise levels below 5 μK-arcmin, this lensing signal becomes an important source of contamination and its sample variance significantly worsens our constraining power on PGWs (Smith et al. 2012). Unlike Galactic foregrounds, the lensing signal is achromatic and cannot be mitigated with multifrequency observations. However, its contribution can be reduced by knowing the cumulative gravitational lensing potential ϕ along the line of sight and having a high-fidelity E-mode map. Together, the two can be combined to form a lensing B-mode template by lensing the E-mode map with the ϕ field and subtracting this template from the measured B-mode map. This technique is known as delensing (Knox & Song 2002; Seljak & Hirata 2004; Carron & Lewis 2017; Carron et al. 2017).

Reconstructing ϕ with high S/N requires high-sensitivity, high-angular-resolution CMB polarization maps (Lewis & Challinor 2006). Therefore, in addition to the low-resolution effort, we assume a separate high-resolution large-aperture instrument dedicated to measuring the intermediate- and small-scale information necessary to delens the B-mode map. This instrument is assumed to have 1' resolution and noise performance equivalent to the 145 GHz channel from the small-aperture telescopes. These experiment specifications are revised for the CMB-S4 Reference Design Report and updated in Section 7. The translation between detector effort and map noise in the delensing instrument is based on the method used for the low-resolution instrument (as described in Section 4.1 and Equation (7)), but using SPTpol achieved performance (Sayre et al. 2020), i.e., without incurring penalties from nonidealities specific to low-resolution instruments and low-ℓ analysis, such as low-ℓ mode filtering and nonuniform coverage. Following the iterative formalism in Smith et al. (2012), ¹⁰¹ using ${{\ell }}_{\min }=300$ and ${{\ell }}_{\max }=4000$ for the ϕ reconstruction and ${{\ell }}_{\min }=30$ for the E-mode map, we convert the map noise in the delensing survey to a delensing efficiency, or equivalently a fractional residual in the lensed B-mode power, specified by setting the residual lensing amplitude A_L to the corresponding level.

The detector effort dedicated to the delensing instrument comes out of the total detector effort budget for the r survey, and the distribution of effort between the low-resolution and delensing instruments is part of the optimization process, as shown in Figures 4 and 5.

Zoom In Zoom Out Reset image size

2.4. Multicomponent Theory Model

2.5. Bandpower Covariance Matrix Rescaling

2.5.1. Signal Scaling

The output model expectation values are also useful in the construction of our bandpower covariance matrix. To construct the BPCM components, we use lensed-ΛCDM + BICEP/Keck noise simulations. However, because we have the individual signal-only, noise-only, and signal × noise terms, we can record all the individual BPCM components:

• sig = signal-only terms Cov(S_i × S_j , S_k × S_l );
• noi = noise-only terms $\mathrm{Cov}({{ \mathcal N }}_{i}\times {{ \mathcal N }}_{j},{{ \mathcal N }}_{k}\times {{ \mathcal N }}_{l}$);
• sn1 = signal × noise terms $\mathrm{Cov}({S}_{i}\times {{ \mathcal N }}_{j},{S}_{k}\times {{ \mathcal N }}_{l}$);
• sn2 = signal × noise terms $\mathrm{Cov}({S}_{i}\times {{ \mathcal N }}_{j},{{ \mathcal N }}_{k}\times {S}_{l}$);
• sn3 = signal × noise terms $\mathrm{Cov}({{ \mathcal N }}_{i}\times {S}_{j},{S}_{k}\times {{ \mathcal N }}_{l}$);
• sn4 = signal × noise terms $\mathrm{Cov}({{ \mathcal N }}_{i}\times {S}_{j},{{ \mathcal N }}_{k}\times {S}_{l}$).

Here, S are signal simulations, ${ \mathcal N }$ are noise simulations, and the indices i, j, k, and l run over the experimental frequency channels.

While calculating the covariances from the signal and noise simulations, we also record the average signal bandpowers from the simulations. For a new signal sky model, we can calculate the new bandpower expectation values, and rescale the signal components in the bandpower covariance matrix by the appropriate power of the ratio of the recorded average signal bandpowers and the newly calculated expectation values. The full BPCM construction and rescaling procedure is presented in Appendix B.

When we do this, we set to zero any term that has an expectation value of zero (under the assumption that signal and noise are uncorrelated, and different signals are uncorrelated) to reduce the Monte Carlo error in the resulting covariance matrix, given the relatively modest number of 499 realizations used. We also set to zero the covariance between bandpowers that are separated by more than one bin in ℓ, but importantly, preserve the covariance between the auto- and cross-spectra of the different frequency bands.

It is worth noting that this procedure allows us to have different numbers of degrees of freedom per bandpower for noise than for signal, which is a byproduct of signal and noise entering differently in a real analysis. This complication is often ignored in other forecasts by setting the noise and signal degrees of freedom to be identical.

The ability to estimate a BPCM for any model means that only a single set of simulations is necessary, and one does not have to run simulations for any and all conceivable scenarios. As already mentioned, in all the projections below we choose to fix Σ(θ) = Σ, and hence we only apply the rescaling step once per fiducial model considered, i.e., we do not rescale our BPCM at every step along the way.

2.5.2. Noise Scaling

In addition to scaling from one signal model to another, recording all the covariance terms allows us to rescale the noise components as well. Given a data set for which we have simulations, the noise scaling can be performed in one of two ways. The first is to take a frequency present in the data set and scale down the noise in the BPCM by the desired amount. In particular, each BPCM component is scaled independently by $\sqrt{{N}_{{\ell },{\rm{S}}4}/{N}_{{\ell },\mathrm{BK}}}$ for each factor of ${ \mathcal N }$ present. The second way is to add an additional frequency, for which we do not yet have data, by taking the covariance structure of an existing frequency, scaling down the noise by the amounts described above, and then expanding the BPCM by filling it in with the appropriate variance and covariance terms between the new band and all the existing ones. These tools facilitate the construction of a new data structure to explore any combination of frequency bands, with any sensitivity in each band.

To obtain N_ℓ values for a CMB-S4 channel by scaling the achieved N_ℓ values of a Stage-3 channel, we have to first scale by the ratio of their respective survey weights and then scale by the ratio of beam window functions, ${B}_{{\ell }}^{2}$. The survey weight is defined as $w=2{f}_{\mathrm{sky}}/{\sigma }_{\mathrm{map}}^{2}$, where f_sky is the effective sky area and σ_map is the Q/U map noise level. The input simulations use measured non-Gaussian B_ℓ shapes, but we rescale based on Gaussian approximations that are close to the true functions. We can write the noise spectrum of a CMB-S4 channel as

$$\begin{eqnarray}&&{N}_{{\ell },{\rm{S}}4}={N}_{{\ell },\mathrm{BK}}\displaystyle \frac{{w}_{\mathrm{BK},\mathrm{achieved}}}{{w}_{{\rm{S}}4,\mathrm{projected}}}\displaystyle \frac{{B}_{{\ell },{\rm{S}}4}^{2}}{{B}_{{\ell },\mathrm{BK}}^{2}},\end{eqnarray} \tag{ 5 }$$

where ${B}_{{\ell },\nu }^{2}=\exp \tfrac{-{\ell }({\ell }+1){{\rm{\Theta }}}_{\nu }^{2}}{8\mathrm{log}(2)}$, Θ_ν is the FWHM, in radians, of the Gaussian beam, and w_i,achieved is the achieved integrated survey weight of a particular instrument.

To obtain the projected achieved survey weight for any of the CMB-S4 channels, we rescale the achieved survey weights as

$$\begin{eqnarray}&&\displaystyle \frac{{w}_{{\rm{S}}4,\mathrm{projected}}}{{w}_{\mathrm{BK},\mathrm{achieved}}}=\displaystyle \frac{{n}_{{\rm{S}}4}^{\det -\mathrm{yr}}}{{n}_{\mathrm{BK}}^{\det -\mathrm{yr}}}\displaystyle \frac{{{NET}}_{\mathrm{BK},\mathrm{ideal}}^{2}}{{{NET}}_{{\rm{S}}4,\mathrm{ideal}}^{2}},\end{eqnarray} \tag{ 6 }$$

where ${n}_{{\rm{S}}4}^{\det -\mathrm{yr}}$ is the number of detector-years assumed for CMB-S4 at any particular frequency, and ${n}_{\mathrm{BK}}^{\det -\mathrm{yr}}$ is the number of detector-years in the BK15 data set, with the instruments in their final state (BICEP2/Keck Array Collaborations X 2018).

The implicit assumption in this step is that the performance of this new CMB-S4 frequency channel falls short of idealized performance by the same factor as the real map from which we are scaling. The survey weight scaling is always performed from the closest frequency for which we have available simulation inputs: 20–95 GHz are scaled from BICEP/Keck 95 GHz data; 145 and 155 GHz are scaled from BICEP/Keck 150 GHz data; and 220 and 270 GHz are scaled from BICEP/Keck 220 GHz data.

Since we are using BICEP/Keck products, which are calculated with a particular sky mask m (with ${f}_{\mathrm{sky}}^{\mathrm{BK}}=\tfrac{{{\rm{\Omega }}}_{\mathrm{pix}}}{4\pi }{\sum }_{i}{m}_{i}\simeq 1$%, where the sum is over the pixel i), we must also scale these products appropriately to evaluate the effect of different sky fractions. We propagate the effects of f_sky in the noise spectra and BPCM in two ways: first, we inflate the N_l values by a factor $\beta ={f}_{\mathrm{sky}}^{{\rm{S}}4}/{f}_{\mathrm{sky}}^{\mathrm{BK}}$, which boosts the (signal × noise) and (noise × noise) terms of the covariance matrix by β and β², to take into account the redistribution of the achieved sensitivity onto a larger patch. Note that the (signal × signal) component remains unchanged in this step. Second, we scale down the entire covariance matrix by a factor of β to increase the number of degrees of freedom in the BPCM, accounting for the fact that we are now observing more modes. This procedure scales the signal and noise degrees of freedom independently, preserving the relative effects that filtering and nonuniform coverage have on the covariance structure. We revisit the way this scaling is performed in Section 7, where we separately take into account the impact of realistic observing strategies on the various components.

With the N_ℓ scalings in hand, we can perform the aforementioned BPCM operations to arrive at a scaled CMB-S4 BPCM that encompasses the intricacies of realistic observing conditions.

In this section, we answer the following question: given a fixed amount of effort and the instrument specifications offered in the previous sections, what is the optimal distribution of effort for foreground cleaning and delensing such that the tightest constraint on r is achieved? To do this, we set up an optimization process that calculates the steepest descent through the ten-dimensional space (effort in the nine single-frequency low-resolution channels plus one high-resolution channel for delensing).

We operate in discrete units of effort, with a single unit defined to be equivalent to 500 detector-years at 150 GHz (similar to a BICEP/Keck 150 GHz receiver observing for 1 yr). For other channels, the number of detector-years per unit of effort is calculated as ${n}_{\det ,150}\times {\left(\nu /150\,\mathrm{GHz}\right)}^{2}$. We define "effort" in these units because it is proportional to focal plane area, which is one of the strongest drivers of the overall project cost.

At each step of the algorithm, we allocate a unit of effort in each dimension. For each separate allocation, we rescale the BPCM, compute a new Fisher matrix, and calculate the resulting σ(r). We then compare the constraints and permanently assign one half of the unit of effort to the channel that produces the largest improvement in σ(r) and the other half to its atmospheric window counterpart (the groupings are 30/40, 85/95, 145/155, and 220/270 GHz). This last step, motivated by earlier iterations through the forecasting loop, is enforcing a split in order to reduce the foreground biases obtained under various foreground models. Projections run to a total of 6000 units of effort, which would be equivalent to 500,000 150 GHz detectors operating for six years. Stage-4 scale surveys seem likely to be in the range of (1–3) × 10⁶ detector-years, an order of magnitude increase from Stage-3 experiments. Though it is generally prohibitive to calculate the entire 10-dimensional hypercube of σ(r), we have validated our approach with a full-grid calculation at various points in the optimization.

Figure 4 shows the optimized constraints on r as a function of total effort, as well as the fraction of effort spent on removing the lensing sample variance and the resulting map rms lensing residual, for the no-detection scenario (i.e., r = 0). ¹⁰² To reach the desired science goal of σ(r) = 5 × 10⁻⁴, for f_sky = 3% (see discussion on sky fraction below), requires 1.2 × 10⁶ 150 GHz equivalent detector-years (or 1.8 × 10⁶ when including marginalization over the decorrelation parameters).

Upon obtaining the optimized detector count distribution, we obtain the input noise spectra according to Equation (5). To use these spectra to create noise simulations (discussed in Section 4.1), it is useful to distill them to a few input parameters. To that end, we fit them to the formula

$$\begin{eqnarray}&&{N}_{{\ell },\mathrm{fit}}=\displaystyle \frac{{\ell }({\ell }+1)}{2\pi }\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{pix}}}{{B}_{{\ell }}^{2}}\left(1+{\left(\displaystyle \frac{{\ell }}{{{\ell }}_{\mathrm{knee}}}\right)}^{\gamma }\right){\sigma }_{\mathrm{map}}^{2}\end{eqnarray} \tag{ 7 }$$

and obtain the map depth σ_map, slope γ, and ℓ_knee values. For the small-aperture data, we find ℓ_knee = 50–60 with γ of −2 to −3, depending on the frequency. The optimal distribution of effort is presented in Figure 5, and the configuration that achieves the science goal in 4 yr (or 6 yr when marginalizing over decorrelation) is summarized in Table 1.

Zoom In Zoom Out Reset image size

As mentioned, it is also necessary to optimize the sky fraction. The trade-off between raw sensitivity, ability to remove foregrounds, and ability to delens results is a complicated optimization problem with respect to sky coverage. Figure 6 shows the r sensitivity forecast for CMB-S4 as a function of the observed sky fraction for the case that we only have an upper limit (r = 0). In this case, the optimization prefers a deep survey that targets as small an area as possible. This conclusion, of course, depends on the forecasting assumptions; to that end, we would like to draw attention to several key factors.

Zoom In Zoom Out Reset image size

First, holding the desired constraint on r fixed, the level to which we rely on delensing to decrease the sample variance increases at smaller sky fractions, as expected. For example, as shown in Figure 6, achieving the forecasted sensitivity on r for a survey targeting 1% of the sky will require a >80% reduction in the map rms level of the CMB lensing B-modes. While from a sensitivity standpoint it is possible to achieve these levels, the extent to which systematic effects and small-scale foregrounds will need to be constrained may become too stringent (Carron et al. 2017; Polarbear Collaboration et al. 2020).

Second, the current optimization assumes uniform foreground behavior across the sky (with amplitude equal to that in the BICEP2/Keck region), while in reality the average amplitude, and possibly the complexity of foregrounds, increase as larger sky areas are targeted. This effect would steepen the optimization curve at larger sky fractions and increase our preference for small amounts of sky.

Third, in the case of a detection, a practical consideration for the robustness of the final r result is its reproducibility across the sky. It is therefore useful to observe multiple 1% patches from which we can derive and compare separate cosmological constraints.

Finally, technical aspects of E/B separation of CMB maps may heavily disfavor patches smaller than about 1% of the sky due to cut-sky effects (Bunn et al. 2003; Smith & Zaldarriaga 2007; BICEP2/Keck Array Collaborations VII 2016).

Balancing the forecasting results with these concerns, we have chosen 3% as the default sky fraction for CMB-S4 r constraints (assuming a true value of r = 0). This choice was made for both the CMB-S4 Science Book (Abazajian et al. 2016) and CMB-S4 CDT Report (Lawrence et al. 2017) versions of the forecasts. As mentioned in Section 2.5.2, we revisit the issue of sky coverage in Section 7 with updated assumptions about our survey strategy and how we perform the f_sky scaling from the achieved products.

Using simulations to optimize the design of a CMB experiment inevitably involves a trade-off between the degree of detail that the simulations are able to capture and the computational cost of generating and analyzing them. This trade-off includes the choice of the domain in which the simulation is generated, ranging from the most detailed (but most expensive) time domain, through the map domain, to the most simplified (but most flexible) spectral domain. Inclusion of additional detail can help validate semianalytic results, explore their sensitivity to assumptions about foreground models, sky coverage, and instrumental noise and systematics, and in more mature stages of design, can inform specific instrument and survey strategy choices.

Here, we review the methods used to explore the parameter space for the PGW survey, including map-level noise simulations, sky models, and observation strategy. We also describe our approach to modeling instrumental systematics, the delensing survey, and the analysis methods.

In addition, we use these simulations to validate the spectral domain forecasts for configurations where the approaches are directly comparable, thereby closing the loop presented in Figure 1.

4.1. Map Noise Simulations

4.2. Foreground Models

4.3. Instrumental Systematics

Control of instrumental systematics is a critical design consideration. However, predicting and modeling these effects realistically is a difficult task that is dependent on actual instrument and survey design details, and furthermore, their impact on actual results comes not through the modeled effects but through unmodeled residuals. Many existing CMB experiments have published in-depth studies that use calibration data and simulations to set upper limits on a wide variety of effects, e.g., Keisler et al. (2015), BICEP2 Collaboration III (2015), Essinger-Hileman et al. (2016), Polarbear Collaboration et al. (2020), and BICEP2/Keck Array Collaborations XI (2019). There have also been several publications that examine the effect of specific classes of instrumental systematics on a generic polarization experiment, e.g., Hu et al. (2003), O'Dea et al. (2007), Shimon et al. (2008), Wallis et al. (2014), and Duivenvoorden et al. (2019).

For this study, in the absence of detailed instrument and survey designs, we have taken the first steps in simulating various generic classes of additive systematic by injecting additional noise-like components into the maps and then reanalyzing them without knowledge of what was put in. We have experimented with components that are both correlated and uncorrelated across frequency bands, and that have white, 1/ℓ, and white +1/ℓ spectra, at varying levels compared to single-frequency map noise, or for correlated cases, combined map noise. The leading-order effects of such components can be mitigated via explicit modeling or filtering, but they may still produce map-level residuals. Examples of mechanisms in this class include bandpass mismatches, beam and pointing variations, calibration variations, cross-talk effects, half-wave plate leakage, ground pickup, and readout irregularities.

To assess the impact of instrumental systematics on measurement requirements, for the purpose of determining both the required survey depths and the maximum allowable levels of systematic effects in the final single-frequency survey maps, our general procedure is to feed parameterizations of various systematic effects into semianalytic forecasts and judge at what levels classes of systematics introduce parameter biases or additional uncertainties that are significant compared to the science targets for those parameters.

For the uncorrelated systematic, we add power independently to each frequency map, scaled by the noise power in that map. This systematic is "noise-like," and therefore the power in the map does not roll off at high multipoles due to beam smoothing. The expression for the auto-spectrum of the systematic is

$$\begin{eqnarray}&&{C}_{l}={\sigma }_{\mathrm{map}}^{2}\left(A+B{\left(\displaystyle \frac{l}{{l}_{\mathrm{knee}}}\right)}^{\gamma }\right),\end{eqnarray} \tag{ 8 }$$

where σ_map, l_knee, and γ are same as found in Section 3; A and B are the amplitudes of the white and 1/ℓ components, in units of fraction of noise power.

For the correlated systematic, we add a common signal at the same amplitude to all frequency maps. Since this signal is common, it affects both auto- and cross-spectra as

$$\begin{eqnarray}&&{C}_{l,{ij}}={\sigma }_{\mathrm{comb}}^{2}\left(A+B{\left(\displaystyle \frac{l}{{l}_{\mathrm{knee}}}\right)}^{\gamma }\right){B}_{l,i}{B}_{l,j}.\end{eqnarray} \tag{ 9 }$$

In place of the per-frequency noise power levels, we define the overall level of this systematic relative to ${\sigma }_{\mathrm{comb}}^{2}$, which is the inverse quadrature sum of the per-frequency noise power levels. We use a common l_knee = 50 and γ = − 2. Since this systematic is "CMB signal-like," we calculate the bias assuming that it has been smoothed in each map according to the beam size of that frequency.

Other classes of systematics could be simulated by manipulating the analysis procedure only. Examples of such effects include uncertainties in the bandpasses, polarization angles, calibration, and beam shapes. Such examples are not presented here and are left for future work.

4.4. Delensing

To make simulated maps, the noise realizations described in Section 4.1 are added to the sky models described in Section 4.2. For each realization, one then has a stack of multifrequency I, Q, and U maps containing nonuniform noise, foregrounds, and signal, and the challenge is to reanalyze them to recover the parameter of interest (in this case, r). This can be done by different teams using different methods, and could be carried out in a blind manner, although we have not done this yet.

So far, we have experimented with two methods. The first is a spectral internal linear combination (ILC) method (e.g., Tegmark et al. 2003), which determines the linear combination of multipole coefficients that minimizes the foreground and noise power without altering the CMB contribution. This method only relies on the frequency dependence of the CMB, and does not rely on assumptions about the spectral dependence of foreground components. Determination of cleaning coefficients in each bin (Δℓ = 31) leads to the smallest foreground residuals, but yields foreground residuals that are difficult to model. Relying on a single cleaning coefficient across all bins leads to foreground residuals that are easy to model but large. As a consequence, we use the same cleaning coefficients in three neighboring bins, which results in residuals that are acceptable and can still be marginalized in the power spectrum likelihood analysis. For the marginalization, since foreground residuals are typically dominated by dust, we assume an ℓ-dependence of the residuals across the three bins that share a common cleaning coefficient that is consistent with the ℓ-dependence observed by Planck for dust. For models that exhibit decorrelation between different dust components, this model is no longer correct. This could be improved with better understanding of foregrounds, but no attempts were made to do so in this analysis method.

The second method is an evolution of the parametric multicomponent fit to the ensemble of auto- and cross-spectra, as used for the BICEP/Keck analysis to date (BICEP2/Keck Array & Planck Collaborations 2015; BICEP2/Keck Array Collaborations VI 2016; BICEP2/Keck Array Collaborations X 2018). This method fits the observed bandpowers to a model composed of the lensing expectation plus dust and synchrotron contributions and a possible r component. Dust and synchrotron each have an amplitude (A_dust and A_sync), a spatial spectral parameter (α_d and α_s), and a frequency spectral parameter (β_d and β_s). We also allow a dust/synchrotron correlation (), and decorrelation of the foreground patterns over frequency (Δ_d and Δ_s). This model is equivalent to the one described in Section 2.4.

Both of these analysis methods are only close to optimal when the foreground behavior is close to uniform across the observing field. For analysis of larger fields, algorithms that fit more complex behavior will likely be required, for example, modeling the frequency spectral indices individually in (large) pixels.

Table 2 summarizes the results of the analysis for simulations of the optimized configuration obtained in Section 3 (1.2 × 10⁶ 150 GHz equivalent detector-years) and residual lensing power A_L = 0.1. ¹⁰⁴ The lensing residual is expected for iterative EB delensing according to Smith et al. (2012) for the sensitivity and angular resolution of the delensing survey. The results from the parametric analysis naturally depend on whether a marginalization over decorrelation is performed, while the ILC analysis did not attempt to capture the effects of decorrelation on the recovery of r and σ(r). This is evidenced by the large bias for the ILC method for Model 5 when compared to the parametric analysis that directly accounts for a possible decorrelation (last column). In general, we see that for r = 0 the simple Gaussian foreground Model 0 gives σ(r) ≈ 5 × 10⁻⁴, exactly as expected from the semianalytic formalism. As we progress to the more complex foreground models, σ(r) is generally in the range (5–8) × 10⁻⁴.

The level of biases is generally below 1.0σ for all the models. These simulations are sets of 499 realizations, so the statistical uncertainty on the bias is approximately 0.04σ. However, the strong decorrelation in Model 5, as well as the high-significance detection of decorrelation in the parametric analysis of Model 4, do significantly increase σ(r) and the level of bias. While the parametric method is able to account for the decorrelation, by construction information is lost, and in fact, if one believed in such a scenario, a different re-optimization to concentrate the sensitivity at closer-in (less decorrelated) frequencies would be called for.

Table 3 shows results on detection significance for the CDT Report configuration for sky Model 6. For r = 0, the 95% upper limit is about 2.1σ(r). The value of the tensor-to-scalar ratio for which we expect a 5σ detection after 4 yr of operation is r = 0.004. ¹⁰⁵ For a tensor-to-scalar ratio of r = 0.003, the median detection significance after 4 yr is expected to be 4σ. If a detection were to be emerging at this point, extending the run time to 8 yr would be justified in order to reach a 5σ detection.

While σ(r) can be precisely forecast for given assumptions, the true achieved detection level for r depends on the particular realization of the B-mode field in the observed patch of sky. Therefore, we can only forecast a distribution of detection levels. For a tensor-to-scalar ratio of r = 0.003 and 8 yr of observing, we expect to achieve more than a 3σ detection with a probability of 0.99, more than 4σ with a probability of 0.93, more than 5σ with a probability of 0.53, and more than 6σ with a probability of 0.14. For simplicity, we focus on σ(r), and on median detection levels as well as median 95% confidence upper limits to state the typical outcome.

The numbers in Table 2 clearly show dependence on the foreground model used in the simulation. If the actual foregrounds are substantially different than any of these cases, then the biases could be larger. To obtain some understanding of how large the biases could be, and what instrument modifications might help to reduce them, we have also looked at ILC biases in the extreme case that the foreground residuals are not modeled or marginalized over, but simply absorbed into the estimated B-mode power spectrum. Doing so with dedicated simulations based on sky Model 6 increases the magnitude of the bias on r to 4.1 × 10⁻⁴. The dominant contribution to the bias comes from synchrotron residuals, which motivated placing a lower-frequency 20 GHz channel on the LAT (see footnote 96), reducing the magnitude of the bias to 1.8 × 10⁻⁴, as shown in Table 2.

Table 4 summarizes the results of the analysis of simulations including additive systematic effects on top of foreground Model 3 (note: these simulations correspond to DC3). Different combinations of uncorrelated and correlated contamination with varying spectra are considered. The levels of systematic contamination for these simulations were chosen to predict biases on r of 1 × 10⁻⁴ in semianalytic forecasts. We can see that the different combinations explored increase the bias on r by amounts that typically vary from 0.5 to 1.5 × 10⁻⁴ for the two separate analyses, over the different cases. We find that, to restrict the bias on r to this level, the sum of additive contamination effects needs to be controlled to 3%−7% of the single-frequency survey noise, or (in the case of correlated systematics) 6%−11% of the total combined noise levels. Such percentages are consistent with the upper limits currently achieved for residual additive systematic contamination compared to survey noise by small-aperture experiments (e.g., BICEP2/Keck Array Collaborations VI 2016). Assuming that CMB-S4 will include a sustained effort to continue to control, understand, and model systematic effects down to levels limited by survey noise, these percentages provide reasonable benchmark requirements.

Results of simulating systematic errors in the determination of bandpasses vary by analysis method. The construction of the ILC method makes it largely insensitive to such errors. The parametric analysis, which includes parameterized models of the frequency spectra of different foregrounds, shows biases on r at the 1 × 10⁻⁴ level for uncorrelated random deviations in band-center determination of 0.8%, or for correlated deviations of 2%; we adopt these as reasonable benchmark requirements to accommodate a variety of both blind and astrophysical foreground modeling approaches.

In this section, we describe updates to the framework and the reference design that attempt to take into account the impact of realistic observing strategies, realistic focal-plane layouts, and mapping onto optics tubes, as well as a more conservative approach to our delensing forecasts. We also seek to answer the question of siting, with the South Pole and Chile being the two choices considered. The work in this section has led to the forecasts and plans presented in the most recent CMB-S4 document—the CMB-S4 Reference Design Report (Abazajian et al. 2019). While this section is self-contained, we direct the reader to the aforementioned source for more details.

In previous versions of our forecasting, we have had the ability to choose the number of detectors in each frequency band in a continuously variable manner, as shown in Figure 5. For the reference design, a mapping of detectors into dichroic optics tubes has been carried out, while seeking to maintain the band distribution as determined in the optimization calculations. In scaling the achieved performance from the existing monochromatic instruments to dichroic detectors and optics, no degradation of optical performance has been assumed at this stage, but this assumption should be verified with upcoming data. This results in the configuration described in Abazajian et al. (2019) with 18 SAT tubes, observing for 7 yr. We use this configuration to scale the BICEP/Keck noise bandpower covariance matrix in the same way as described in Section 2.5. A further rescaling is applied to account for sky coverage, as explained below.

7.1. Sky Coverage Effects

The semianalytic calculations of Section 2.5 assumed a simplified rescaling for sky area, while the map-based simulations of Section 6 assumed an idealized circular sky patch, which is not actually achievable with a practical instrument from a site at any latitude. Figure 7 compares our prior assumptions to more realistic hit patterns.

Zoom In Zoom Out Reset image size

From the South Pole, it is possible to concentrate the coverage onto a compact region of sky, but from Chile, the region that can be observed is affected by Earth's rotation, resulting in more extended coverage. The large instantaneous field of view of the SAT telescopes means that there is a minimum field size that can be achieved, and that there is always a strong edge taper in the coverage pattern.

We have performed a calculation that attempts to optimize simulated SAT observations from Chile to produce the densest possible coverage on a 3% patch of low-foreground sky, resulting in the overall pattern shown in Figure 7 as "Chile full." We segment this into its deepest part, which we call "Chile deep," and the remainder, which we call "Chile shallow."

From the South Pole, one can scan the same patch at all times of the day and year at the same observing elevation, with the size of the observed patch controlled by the length of the scan throw in R.A. A minimal-length scan results in the pattern shown in the figure as "Pole deep." Lengthening the scan while remaining in the low foreground sky results in the pattern "Pole wide." In the results below, "Pole deep" and "Pole wide" are therefore "either/or" options.

Because the noise increases in regions with less observing time, the effective sky area for noise is larger than the effective sky area for signal—and both of these also depend on the weighting applied when analyzing the maps. The patterns shown in Figure 7 have the effective sky fractions reported in Table 5, assuming inverse noise-variance weighting.

Table 5. Effective Sky Fractions for Signal, Noise, and Signal-cross-noise, as Percentages, for the Observation Patterns Shown in Figure 7, and the Case of Inverse Noise-variance Weighting (i.e., Equations (11)–(13) with w_i = h_i )

	SP Deep	SP Wide	CH Deep	CH Shallow	CH Full
${f}_{\mathrm{sky}}^{\mathrm{sig}}$	1.9	4.3	2.4	10	5.9
${f}_{\mathrm{sky}}^{\mathrm{noi}}$	2.9	6.5	3.4	20	18
${f}_{\mathrm{sky}}^{\mathrm{cross}}$	2.5	5.5	3.0	16	12

Notes. Here, "SP" is the South Pole and "CH" is Chile. The CH Shallow numbers appear larger than CH Full due to the effects of the weighting.

For comparison, the BICEP2/Keck and BICEP3 values are (1.0, 1.3, and 1.1) and (1.9, 2.7, and 2.3) for signal, noise, and signal-cross-noise respectively.

Download table as:  ASCIITypeset image

We can take into account the above effects by rescaling the BICEP/Keck BPCM in a more sophisticated manner. First, we need to scale the noise due to distributing the effort on a patch of sky larger than the original BICEP/Keck one. The noise is scaled by the effective noise factor

$$\begin{eqnarray}&&{f}_{\mathrm{eff}}^{\mathrm{noise}}=\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{pix}}}{4\pi }\displaystyle \frac{{\sum }_{i}{w}_{i}^{2}{h}_{i}^{-1}{\sum }_{i}{h}_{i}}{{\sum }_{i}{w}_{i}^{2}},\end{eqnarray} \tag{ 10 }$$

where Ω_pix is the solid angle of a single pixel, w_i are the weights for pixel i, and h_i are the hit counts.

Second, since we observe a different number of modes, we need to scale the signal, noise, and signal-cross-noise contributions of the BPCM by the factors

$$\begin{eqnarray}&&{f}_{\mathrm{sky}}^{\mathrm{noise}}=\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{pix}}}{4\pi }\displaystyle \frac{{\left({\sum }_{i}{w}_{i}^{2}{h}_{i}^{-1}\right)}^{2}}{{\sum }_{i}{w}_{i}^{4}{h}_{i}^{-2}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{f}_{\mathrm{sky}}^{\mathrm{signal}}=\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{pix}}}{4\pi }\displaystyle \frac{{\left({\sum }_{i}{w}_{i}^{2}\right)}^{2}}{{\sum }_{i}{w}_{i}^{4}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{f}_{\mathrm{sky}}^{\mathrm{cross}}=\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{pix}}}{4\pi }\displaystyle \frac{{\sum }_{i}{w}_{i}^{2}\sum _{i}{w}_{i}^{2}{h}_{i}^{-1}}{{\sum }_{i}{w}_{i}^{4}{h}_{i}^{-1}}.\end{eqnarray} \tag{ 13 }$$

We also need to take out the effect of these factors from the original BPCM. In the BICEP/Keck analysis, the weights are the inverse noise variance, i.e., w_i = h_i . In the case of CMB-S4, we will never be noise-dominated, either due to an actual primordial signal, or due to the lensing and foreground residuals, so here we use the inverse variance of the total signal and noise to determine the optimal weighting.

Last, while the scanning strategy used here has been optimized to mostly observe the cleanest available sky, some regions with high Galactic emission are also observed. Realistically, we would mask such regions when analyzing the maps. To assess this effect, we used masks based on a smoothed Planck 353 GHz polarized intensity map, keeping the cleanest 30% or 60% of the full sky (28% and 58% after apodization). We then used these masks to disregard some of the pixels, resulting in a degradation of the constraints on r.

7.2. Delensing Revised

7.3. Results

The covariance matrices calculated as described above are used to produce the results given in this section, where the number, siting, and coverage patterns of the SATs are varied. In all cases, a delensing LAT at the South Pole is assumed to concentrate its coverage on a small patch of sky, while delensing over larger sky areas is assumed to be available from the Chilean LATs.

As mentioned earlier, we split the Chilean coverage shown as "Chile full" in Figure 7 into a deep patch, "Chile deep," which overlaps with the "Pole deep" region, and call the remainder "Chile shallow." We then make separate forecasts for each Chilean subregion using the appropriate delensing level for that subregion. To obtain forecasts for the full region, we add the σ(r) results in simple inverse quadrature, thereby making the approximation of independence of the measured modes. When we combine these with South Pole observations, we mimic a joint analysis over the overlapping region by taking the sum of the "Pole deep" and "Chile deep" coverage maps and computing the corresponding weights and lensing residuals, and add the "Chile shallow" results in inverse quadrature.

The foreground fiducial model is kept the same as in Section 2.4, and we also explore the possibility of unmodeled foreground residuals contributing residual power to the cleaned maps, assuming a foreground cleaning to 1% at 95 GHz. We do so by adding in quadrature to the r forecasts a foreground bias equal to 1% of the equivalent r of the foreground minimum of the BK15 data at 95 GHz and ℓ = 80.

Since some parts of the "Chile shallow" coverage lie closer to the Galactic plane, we boost the foreground level in the shallow region, and the foreground bias if it is included, by a factor of 3 with respect to the deep patch. This scaling is based on the dust amplitude measured in the map-based simulations for surveys observed from the South Pole or from Chile, using a map with spatially varying foregrounds based on Planck data.

In Figure 8, we show the dependence of σ(r) on r for the different coverage masks. We calculate these constraints for r = 0, 0.003, 0.01, and 0.03, for the different Galactic cuts, and show the linear interpolation between these points, sampled on a high-resolution linearly spaced grid. We find that the survey strategy from the South Pole is always favored in the limit of small r, with the crossover point depending on the specific assumptions. This comes as a direct result of the fact that, while a larger fraction of the sky can be observed from Chile, one can concentrate the available sensitivity more deeply from the South Pole (see Figure 7). We note that the delensing requirements are more stringent for the smaller surveys (such as the South Pole one), highlighting the importance of periodically revisiting the assumptions made for the delensing survey of CMB-S4 as new data from Stage-3 experiments become available.

Zoom In Zoom Out Reset image size

In Tables 6–7, we present a set of σ(r) results for seven years of observations and r = 0, and in Tables 8–9, we present a set of significance of detection levels for r = 0.003, while varying the number of optics tubes at the South Pole and Chile over a wide range of options. We show results for two different variants: (i) with no marginalization over the decorrelation parameters and for the 28% cleanest polarized sky; and (ii) the same as (i), but marginalizing over the foreground decorrelation parameters.

Table 6. Combined 10⁴ × σ(r) Values (Smaller Numbers are Better), Assuming r = 0 after 7 yr of Observation, Keeping only the 28% Cleanest Part of the Sky, Assuming no Decorrelation and an Observing Efficiency in Chile the Same as at the South Pole

Note. The bolded cells conserve the 18 SAT optics tube count of the reference design, while nonbolded cells explore other counts at the South Pole and in Chile.

Download table as:  ASCIITypeset image

Table 7. Same as Table 6, but Assuming Additional Foreground Decorrelation Parameters

Download table as:  ASCIITypeset image

Table 8. Combined Detection Significance (Larger Numbers are Better) for r = 0.003 after 7 yr of Observation, Keeping only the 28% Cleanest Part of the Sky, Assuming no Decorrelation and an Observing Efficiency in Chile the Same as at the South Pole

Note. The bolded cells conserve the 18 SAT optics tube count of the reference design, while nonbolded cells explore other counts at the South Pole and in Chile.

Download table as:  ASCIITypeset image

In all these calculations, we assume that the observing efficiency from Chile is 100% equivalent to that from the South Pole. Despite long histories of CMB observations at both sites, it is still quite challenging to make a clean comparison of their observing efficiencies. Results from the BICEP/Keck program are responsible for the leading constraints on r for the last decade, but it is still not currently possible to disentangle the role of the observing site from other factors that have contributed to the success of that program, such as detector performance, instrument design, observing strategy, and operations management. However, note also that even assuming equal observing efficiency at both sites (Table 6), there is still a difference between σ(r) obtained for equivalent numbers of optics tubes in Chile versus the South Pole, due to sky fraction, foreground complexity, and delensing. In the CMB-S4 Reference Design Report, we also show forecasting results that assume an observing efficiency from Chile equivalent to 50% of the South Pole efficiency.

Note that the results where the total number of SAT optics tubes is different from the assumed number of 18 are subject to a caveat: the delensing effort is assumed to be held fixed to the reference design—one LAT at the South Pole and two in Chile. In principle, as the total effort is varied away from the reference design, one should re-optimize the fraction of delensing effort as per Figure 4.

In this paper, we have presented the CMB-S4 r forecasting framework and its evolution through three major iterations. To ensure forecasting realism, we have chosen to directly ground our forecasting in the end-to-end on-sky achieved performance of Stage 2–3 experiments, rather than rely on ab initio assumptions. Additionally, we have checked the semianalytic calculations with map-based simulations of increasing complexity, thereby creating a closed forecasting loop (presented in Figure 1). This approach allows us to confidently incorporate all the effects that impact current real-world experiments, to flexibly trace the end-to-end effects of changes to experimental design or assumptions, and to iterate over the framework as more sophisticated sky and instrument models become available.

The semianalytic calculations described in Sections 2 and 3 indicate that, for a 3% sky fraction, 1.2 × 10⁶ 150 GHz equivalent detector-years (or 1.8 × 10⁶ when including marginalization over the decorrelation parameters) are required to reach the science requirement of σ(r) = 5 × 10⁻⁴, with 30% of this effort assigned to the delensing survey. The resulting optimal distribution of detectors across frequency bands is summarized in Table 3.

The map-based simulations described in Sections 4–6 confirm the σ(r) results from the semianalytic calculations. These simulations also indicate that bias in the recovered r value is within 1σ for a suite of different foreground models. However, we note that foregrounds remain a serious issue that must be periodically revisited as the project progresses.

Additional map-based simulations indicate that the systematic bias on r can be controlled to <1σ, provided that fractional contamination levels similar to those already achieved by small-aperture telescopes can be maintained.

Mapping the requirements from the semianalytic calculations onto realizable instruments (see Section 7) results in the reference design described in Abazajian et al. (2019) for a seven-year survey period, and in the constraints presented in Tables 6–9.

Table 9. Same as Table 8, but Assuming Additional Foreground Decorrelation Parameters

Download table as:  ASCIITypeset image

In Section 7, we have also extended the semianalytic calculations to account for realistic observation patterns and probed the dependence of σ(r) on r, as a function of experiment siting, finding that the survey strategy from the South Pole is always favored in the limit of small r. We note that surveys with smaller footprints (such as the South Pole survey) depend more tightly on the levels of achieved delensing. Therefore, revising the delensing assumptions will be important as new studies of small-scale foregrounds become available.

Going forward, as the CMB-S4 project matures, the collaboration will need to converge on increasingly specific instrumentation choices. To quantitatively discern between the different options, we intend to continue using and developing the closed-loop framework presented here. Among others, we anticipate adding complexity to the following directions: delensing treatment; foreground simulations; foreground parameterization; survey strategies; and instrumental systematics. Additionally, we expect to apply new analysis methods to our data challenges and incorporate new achieved performance levels at multiple additional frequencies from multiple sites as these data become available. The iterative nature of our framework can easily accommodate these revisions, achieving increased realism with each iteration.

The CMB-S4 collaboration (https://cmb-s4.org/) is working to plan, construct, and operate a next-generation, multisite CMB experiment in the 2020s. The collaboration is led by an elected Governing Board, Spokespeople, Committee Chairs, and Executive Team. Funding for the CMB-S4 Integrated Project Office is provided by the Department of Energy's Office of Science (project level CD-0) and by the National Science Foundation through the Mid-Scale Research Infrastructure-R1 award OPP-1935892. This research used resources of Argonne National Laboratory, a U.S. Department of Energy (DOE) Office of Science User Facility operated under Contract No. DE-AC02-06CH11357. This document was prepared by the CMB-S4 collaboration using the resources of the Fermi National Accelerator Laboratory (Fermilab), a U.S. Department of Energy, Office of Science, HEP User Facility. Fermilab is managed by Fermi Research Alliance, LLC (FRA), acting under Contract No. DE-AC02-07CH11359. Work at Lawrence Berkeley National Laboratory was supported by the Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Work at SLAC National Accelerator Laboratory is supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Contract No. DE-AC02-76SF00515. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. In the United States, work on CMB-S4 by individual investigators has been supported by the National Science Foundation (awards 1248097, 1255358, 1815887, 1835865, 1852617, 2009469), the Department of Energy (awards DE-SC0009919, DE-SC0009946, DE-SC0010129, DE-SC0011784), and the National Aeronautics and Space Administration (award ATP-80NSSC20K0518). In Australia, the Melbourne authors acknowledge support from an Australian Research Council Future Fellowship (FT150100074). In Canada, R.H. is supported by the Discovery Grants program from NSERC, and acknowledges funding from CIFAR, the Sloan Foundation, and the Dunlap family. In Italy, C.B. acknowledges support under the ASI COSMOS and INFN INDARK programs. In the Netherlands, D.M. acknowledges NWO VIDI award number 639.042.730. In Switzerland, J.C. is supported by an SNSF Eccellenza Professorial Fellowship (No. 186879). In the United Kingdom, A.L., G.F., and J.C. are supported by the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC grant Agreement No. [616170]. A.L. also acknowledges STFC award ST/P000525/1. S.M. is supported by the research program Innovational Research Incentives Scheme (Vernieuwingsimpuls), which is financed by the Netherlands Organization for Scientific Research through the NWO VIDI grant No. 639.042.612-Nissanke and the Labex ILP (reference ANR-10-LABX-63) part of the Idex SUPER, received financial state aid managed by the Agence Nationale de la Recherche, as part of the program Investissements d'avenir under the reference ANR-11-IDEX-0004-02. Some computations in this paper were run on the Odyssey cluster, supported by the FAS Science Division Research Computing Group at Harvard University.

The forecasting framework uses a parametric model to describe the bandpower expectation values as a combination of cosmological and foreground signals. This parameterization follows the one presented in Appendix G of BICEP2/Keck Array Collaborations X (2018).

The model includes signals from the CMB (lensed-scalar and tensor contributions), Galactic dust, and Galactic synchrotron, with the possibility of spatial correlation between dust and synchrotron. Contributions from dust and synchrotron to the BB spectrum between maps at frequencies ν₁ and ν₂ can be written in the following form, after computing the expected spectrum as a function of ℓ, and bandpower window functions are applied to calculate binned bandpower expectation values:

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal D }}_{{\ell },{BB}}^{{\nu }_{1}\times {\nu }_{2}} & = & {A}_{\mathrm{dust}}{{\rm{\Delta }}}_{{\rm{d}}}^{{\prime} }{f}_{{\rm{d}}}^{{\nu }_{1}}{f}_{{\rm{d}}}^{{\nu }_{2}}{\left(\displaystyle \frac{{\ell }}{80}\right)}^{{\alpha }_{{\rm{d}}}}\\ & & +\,{A}_{\mathrm{sync}}{\rm{\Delta }}{{\prime} }_{{\rm{s}}}{f}_{{\rm{s}}}^{{\nu }_{1}}{f}_{{\rm{s}}}^{{\nu }_{2}}{\left(\displaystyle \frac{{\ell }}{80}\right)}^{{\alpha }_{{\rm{s}}}}\\ & & +\,\epsilon \sqrt{{A}_{\mathrm{dust}}{A}_{\mathrm{sync}}}({f}_{{\rm{d}}}^{{\nu }_{1}}{f}_{{\rm{s}}}^{{\nu }_{2}}+{f}_{{\rm{s}}}^{{\nu }_{1}}{f}_{{\rm{d}}}^{{\nu }_{2}}){\left(\displaystyle \frac{{\ell }}{80}\right)}^{({\alpha }_{{\rm{d}}}+{\alpha }_{{\rm{s}}})/2}.\end{array}\end{eqnarray} \tag{ A1 }$$

Parameters A_dust and A_sync specify the dust and synchrotron power in units of $\mu {{\rm{K}}}_{\mathrm{CMB}}^{2}$ at angular scale ℓ = 80. These are defined at pivot frequencies of 353 GHz for dust and 23 GHz for synchrotron. The dust and synchrotron components scale as power laws in ℓ with slopes α_d and α_s, respectively (note that slope parameters are defined for ${{ \mathcal D }}_{{\ell }}\equiv {\ell }({\ell }+1){C}_{{\ell }}/2\pi$). The parameter specifies the level of spatial correlation between dust and synchrotron; this correlation coefficient is assumed to be constant across all ℓ. If either A_dust or A_sync are negative, the contribution of the correlated component to the expectation value flips sign. Negative foreground amplitudes are technically nonphysical, but this analytic continuation becomes important when we explore entire parameter phase spaces—for instance, in the simulation reanalysis steps.

Coefficients ${f}_{{\rm{d}}}^{\nu }$ and ${f}_{{\rm{s}}}^{\nu }$, given by Equations (A5) and (A6), describe the scaling of dust and synchrotron amplitude from the pivot frequencies to the actual bandpasses of the maps at frequencies ν₁ and ν₂. The SED model used for dust emission is a modified blackbody, i.e., a Planck function multiplied by a power law with emissivity spectral index β_d (Planck Collaboration XXII 2015). We adopt a dust temperature T_d = 19.6 K. The SED model used for synchrotron is defined as a power law with spectral index β_s relative to a Rayleigh–Jeans spectrum. The frequency scaling also includes the unit conversion for μK_CMB from the pivot frequency to the target bandpass.

To integrate the SED and unit conversion over the bandpass of the target map, we adopt the convention used by Planck Collaboration IX (2014), in which our bandpass functions describe response as a function of frequency to a beam-filling source with uniform spectral radiance. For emissivity spectral index β_d and dust temperature T_d, the modified blackbody scaling (MBBS) from pivot frequency ν_pivot to a map with bandpass R(ν) is given by

$$\begin{eqnarray}&&{f}_{\mathrm{MBBS}}^{\nu }=\displaystyle \frac{\int d\nu R(\nu ){\nu }^{3+{\beta }_{{\rm{d}}}}{\left(\exp \tfrac{h\nu }{{{kT}}_{{\rm{d}}}}-1\right)}^{-1}}{{\nu }_{\mathrm{pivot}}^{3+{\beta }_{{\rm{d}}}}{\left(\exp \tfrac{h{\nu }_{\mathrm{pivot}}}{{{kT}}_{{\rm{d}}}}-1\right)}^{-1}}.\end{eqnarray} \tag{ A2 }$$

For a synchrotron power-law scaling (PLS) with spectral index β_s, we calculate the coefficient in a similar way:

$$\begin{eqnarray}&&{f}_{\mathrm{PLS}}^{\nu }=\displaystyle \frac{\int d\nu R(\nu ){\nu }^{2+{\beta }_{{\rm{s}}}}}{{\nu }_{\mathrm{pivot}}^{2+{\beta }_{{\rm{s}}}}}.\end{eqnarray} \tag{ A3 }$$

The conversion between μK_CMB units at the pivot frequency and μK_CMB units at the target map bandpass is given by

$$\begin{eqnarray}&&{f}_{{{\rm{T}}}_{\mathrm{CMB}}}^{\nu }=\displaystyle \frac{{\nu }_{\mathrm{pivot}}^{4}\exp \tfrac{h{\nu }_{\mathrm{pivot}}}{{{kT}}_{\mathrm{CMB}}}{\left(\exp \tfrac{h{\nu }_{\mathrm{pivot}}}{{{kT}}_{\mathrm{CMB}}}-1\right)}^{-2}}{\int d\nu R(\nu ){\nu }^{4}\exp \tfrac{h\nu }{{{kT}}_{\mathrm{CMB}}}{\left(\exp \tfrac{h\nu }{{{kT}}_{\mathrm{CMB}}}-1\right)}^{-2}}.\end{eqnarray} \tag{ A4 }$$

Combining these factors, we obtain the scalings used in Equation (A1):

$$\begin{eqnarray}&&{f}_{{\rm{d}}}^{\nu }={f}_{\mathrm{MBBS}}^{\nu }\times {f}_{{T}_{\mathrm{CMB}}}^{\nu },\end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray}&&{f}_{{\rm{s}}}^{\nu }={f}_{\mathrm{PLS}}^{\nu }\times {f}_{{T}_{\mathrm{CMB}}}^{\nu }.\end{eqnarray} \tag{ A6 }$$

We also consider dust and synchrotron frequency decorrelation. The simplest possible model of a polarized foreground component is one with a fixed spatial pattern on the sky that scales with frequency according to a single SED. In this case, the expectation value of the cross-spectrum between any two frequencies is the geometric mean of the respective auto-spectrum expectation values. In reality, the spatial pattern might vary as a function of frequency, leading to the cross-spectra being suppressed with respect to the geometric mean of the auto-spectra (Planck Collaboration 2017). We refer to this phenomenon as decorrelation.

We model decorrelation in the power spectrum domain using a set of simple one-parameter models. We define the correlation ratio of dust between two reference frequencies, 217 and 353 GHz, at pivot scale ℓ = 80 as

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{d}}}=\displaystyle \frac{{{ \mathcal D }}_{80}(217\times 353)}{\sqrt{{{ \mathcal D }}_{80}(217\times 217){{ \mathcal D }}_{80}(353\times 353)}},\end{eqnarray} \tag{ A7 }$$

where ${{ \mathcal D }}_{80}$ is the dust power at ℓ = 80. Here, Δ_d < 1 corresponds to decorrelation. We scale to other frequency combinations using the factor suggested by Planck Collaboration (2017):

$$\begin{eqnarray}&&f({\nu }_{1},{\nu }_{2})=\displaystyle \frac{{\left(\mathrm{log}({\nu }_{1}/{\nu }_{2})\right)}^{2}}{{\left(\mathrm{log}(217/353)\right)}^{2}}.\end{eqnarray} \tag{ A8 }$$

Similarly, based on suggestions from Planck Collaboration (2017), we consider three possible scalings with ℓ:

$$\begin{eqnarray}g({\ell })=\left\{\begin{array}{ll}1, & \mathrm{flat}\ \mathrm{case};\\ ({\ell }/80), & \mathrm{linear}\ \mathrm{case};\\ {\left({\ell }/80\right)}^{2}, & \mathrm{quadratic}\ \mathrm{case}.\end{array}\right.\end{eqnarray} \tag{ A9 }$$

The scalings above can produce extreme (and nonphysical) behavior at high ℓ or for frequencies that are widely separated. We therefore remap the scaled correlation coefficient using the following function:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{d}}}^{{\prime} }({\nu }_{1},{\nu }_{2},{\ell })=\exp \left[\mathrm{log}({{\rm{\Delta }}}_{{\rm{d}}})\,f({\nu }_{1},{\nu }_{2})\,g({\ell })\right].\end{eqnarray} \tag{ A10 }$$

With this remapping, ${{\rm{\Delta }}}_{{\rm{d}}}^{{\prime} }$ remains in the range 0 (no correlation) to 1 (perfect correlation) for all values of f and g. This combination of frequency scaling and nonlinear remapping has been shown to correspond to a Gaussian spatial variation in the foreground spectral index parameter (Vansyngel et al. 2017).

In a similar vein, we define the parameter ${\rm{\Delta }}{{\prime} }_{{\rm{s}}}$, which describes decorrelation of the synchrotron pattern. We do not include foreground decorrelation parameters in the dust–synchrotron correlated component. A complete foreground model would include all correlations between dust and synchrotron foregrounds across observing frequencies, but the current data do not provide useful guidance about the form of such correlations.

In addition to foregrounds, we include CMB scalar and tensor contributions. We make the simplifying assumptions that the tensor BB spectrum is given by a template scaled by parameter r and the BB spectrum from lensed CMB is given by a template scaled by parameter A_L. Using CMB temperature units, the CMB contribution to the BB spectrum is given by

$$\begin{eqnarray}&&{D}_{{\ell },{BB}}^{{\nu }_{1}\times {\nu }_{2}}=\displaystyle \frac{r}{0.1}{D}_{{\ell },{BB}}^{\mathrm{tensor}}+{A}_{{\rm{L}}}{D}_{{\ell },{BB}}^{\mathrm{lensing}},\end{eqnarray} \tag{ A11 }$$

where ${D}_{{\ell },{BB}}^{\mathrm{tensor}}$ is the BB template for a tensor signal with r = 0.1 and ${D}_{{\ell },{BB}}^{\mathrm{lensing}}$ is the expected lensing BB spectrum for ΛCDM. These are obtained using the CAMB ¹⁰⁸ package (Lewis & Challinor 2011).

The BPCM construction relies on analytic rescaling of simulations. We calculate and store the bandpower covariance of signal and noise simulations for a particular sky and instrument model; rescaling from these covariance matrices eliminates the computational cost of creating large suites of simulations for each desired model. In the rescaling process, we explicitly set to zero terms that have zero expectation value (under the assumption that signal and noise are uncorrelated, and that physically unrelated signals are also uncorrelated) to reduce Monte Carlo error in the covariance matrix that stems from the modest number (499) of simulation realizations used. This framework is described in detail in Buza (2019), but we review it briefly here.

In a general case, let m_a , m_b , m_c , and m_d denote the four maps included in our analysis, with m_a × m_b denoting the cross-spectrum between those maps. Each map is the sum of independent signal components (CMB and foregrounds) plus a noise contribution, i.e., m_a = ∑_i s_ai + n_a . Since the synchrotron and dust foregrounds could be spatially correlated, we choose to divide them into three mutually independent components—uncorrelated synchrotron, uncorrelated dust, and the correlated part of synchrotron and dust. Using the properties that these signal and noise fields are independent and have zero mean, the expectation values for many spectra can be set to zero: $\left\langle {s}_{{xi}}\times {s}_{{yj}}\right\rangle =0$ for i ≠ j (independent signal fields) and $\left\langle {s}_{{xi}}\times {n}_{y}\right\rangle =0$ (signal is independent of noise). Then, a generic bandpower covariance term can be written as

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{Cov}\left({m}_{a}\times {m}_{b},{m}_{c}\times {m}_{d}\right)=\left\langle ({m}_{a}\times {m}_{b})({m}_{c}\times {m}_{d})\right\rangle \\ \ \ \ \ \ -\,\left\langle {m}_{a}\times {m}_{b}\right\rangle \left\langle {m}_{c}\times {m}_{d}\right\rangle \end{array}\end{eqnarray} \tag{ B1 }$$

$$\begin{eqnarray}\begin{array}{rcl} & = & \sum _{i}\mathrm{Cov}\left({s}_{{ai}}\times {s}_{{bi}},{s}_{{ci}}\times {s}_{{di}}\right)\\ & & +\,\sum _{i}\sum _{j\ne i}\mathrm{Cov}\left({s}_{{ai}}\times {s}_{{bj}},{s}_{{ci}}\times {s}_{{dj}}\right)\\ & & +\,\sum _{i}\sum _{j\ne i}\mathrm{Cov}\left({s}_{{ai}}\times {s}_{{bj}},{s}_{{cj}}\times {s}_{{dj}}\right)\\ & & +\,\sum _{i}\mathrm{Cov}\left({s}_{{ai}}\times {n}_{b},{s}_{{ci}}\times {n}_{d}\right)\\ & & +\,\sum _{i}\mathrm{Cov}\left({s}_{{ai}}\times {n}_{b},{n}_{c}\times {s}_{{di}}\right)\\ & & +\,\sum _{i}\mathrm{Cov}\left({n}_{a}\times {s}_{{bi}},{s}_{{ci}}\times {n}_{d}\right)\\ & & +\,\sum _{i}\mathrm{Cov}\left({n}_{a}\times {s}_{{bi}},{n}_{c}\times {s}_{{di}}\right)\\ & & +\,\mathrm{Cov}\left({n}_{a}\times {n}_{b},{n}_{c}\times {n}_{d}\right),\end{array}\end{eqnarray} \tag{ B2 }$$

where the first three terms are signal-only covariances, the next four terms are covariances between select signal-noise cross-spectra, and the last term is a noise-only covariance. Additional terms can be set to zero if we make the further assumption that the noise is independent in each map, i.e., $\left\langle {n}_{x}\times {n}_{y}\right\rangle =0$ for x ≠ y. We also set to zero the covariance between any bandpowers that are separated by more than one bin in ℓ (for bins with Δℓ = 35), since these correlations are very small and not well-measured.

Rather than running simulations with a complicated combination of CMB and foreground signals, we calculate the above covariance terms for a simple signal model and then rescale to obtain the bandpower covariance matrix for an arbitrary model. In practice, the signal simulations are CMB realizations for a lensed ΛCDM model. Using s_x0 to denote the simulated signal component of map m_x , the rescaled bandpower covariance matrix term is given by

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{Cov}\left({m}_{a}\times {m}_{b},{m}_{c}\times {m}_{d}\right)\\ \ \ =\ \mathrm{Cov}\left({s}_{a0}\times {s}_{b0},{s}_{c0}\times {s}_{d0}\right)\\ \ \ \times \,\sum _{i}{\left[\displaystyle \frac{\left\langle {s}_{{ai}}\times {s}_{{ai}}\right\rangle \left\langle {s}_{{bi}}\times {s}_{{bi}}\right\rangle \left\langle {s}_{{ci}}\times {s}_{{ci}}\right\rangle \left\langle {s}_{{di}}\times {s}_{{di}}\right\rangle }{\left\langle {s}_{a0}\times {s}_{a0}\right\rangle \left\langle {s}_{b0}\times {s}_{b0}\right\rangle \left\langle {s}_{c0}\times {s}_{c0}\right\rangle \left\langle {s}_{d0}\times {s}_{d0}\right\rangle }\right]}^{1/2}\\ \ \ +\,\displaystyle \frac{1}{2}\mathrm{Cov}\left({s}_{a0}\times {s}_{b0},{s}_{c0}\times {s}_{d0}\right)\\ \ \ \times \,\sum _{i}\sum _{j\ne i}{\left[\displaystyle \frac{\left\langle {s}_{{ai}}\times {s}_{{ai}}\right\rangle \left\langle {s}_{{bj}}\times {s}_{{bj}}\right\rangle \left\langle {s}_{{ci}}\times {s}_{{ci}}\right\rangle \left\langle {s}_{{dj}}\times {s}_{{dj}}\right\rangle }{\left\langle {s}_{a0}\times {s}_{a0}\right\rangle \left\langle {s}_{b0}\times {s}_{b0}\right\rangle \left\langle {s}_{c0}\times {s}_{c0}\right\rangle \left\langle {s}_{d0}\times {s}_{d0}\right\rangle }\right]}^{1/2}\\ \ \ +\,\displaystyle \frac{1}{2}\mathrm{Cov}\left({s}_{a0}\times {s}_{b0},{s}_{c0}\times {s}_{d0}\right)\\ \ \ \times \,\sum _{i}\sum _{j\ne i}{\left[\displaystyle \frac{\left\langle {s}_{{ai}}\times {s}_{{ai}}\right\rangle \left\langle {s}_{{bj}}\times {s}_{{bj}}\right\rangle \left\langle {s}_{{cj}}\times {s}_{{cj}}\right\rangle \left\langle {s}_{{di}}\times {s}_{{di}}\right\rangle }{\left\langle {s}_{a0}\times {s}_{a0}\right\rangle \left\langle {s}_{b0}\times {s}_{b0}\right\rangle \left\langle {s}_{c0}\times {s}_{c0}\right\rangle \left\langle {s}_{d0}\times {s}_{d0}\right\rangle }\right]}^{1/2}\\ \ \ +\,\mathrm{Cov}\left({s}_{a0}\times {n}_{b},{s}_{c0}\times {n}_{d}\right)\\ \ \ \times \,\sum _{i}{\left[\displaystyle \frac{\left\langle {s}_{{ai}}\times {s}_{{ai}}\right\rangle \left\langle {s}_{{ci}}\times {s}_{{ci}}\right\rangle }{\left\langle {s}_{a0}\times {s}_{a0}\right\rangle \left\langle {s}_{c0}\times {s}_{c0}\right\rangle }\right]}^{1/2}\\ \ \ +\,\mathrm{Cov}\left({s}_{a0}\times {n}_{b},{n}_{c}\times {s}_{d0}\right)\\ \ \ \times \,\sum _{i}{\left[\displaystyle \frac{\left\langle {s}_{{ai}}\times {s}_{{ai}}\right\rangle \left\langle {s}_{{di}}\times {s}_{{di}}\right\rangle }{\left\langle {s}_{a0}\times {s}_{a0}\right\rangle \left\langle {s}_{d0}\times {s}_{d0}\right\rangle }\right]}^{1/2}\\ \ \ +\,\mathrm{Cov}\left({n}_{a}\times {s}_{b0},{s}_{c0}\times {n}_{d}\right)\\ \ \ \times \,\sum _{i}{\left[\displaystyle \frac{\left\langle {s}_{{bi}}\times {s}_{{bi}}\right\rangle \left\langle {s}_{{ci}}\times {s}_{{ci}}\right\rangle }{\left\langle {s}_{b0}\times {s}_{b0}\right\rangle \left\langle {s}_{c0}\times {s}_{c0}\right\rangle }\right]}^{1/2}\\ \ \ +\,\mathrm{Cov}\left({n}_{a}\times {s}_{b0},{n}_{c}\times {s}_{d0}\right)\\ \ \ \times \,\sum _{i}{\left[\displaystyle \frac{\left\langle {s}_{{bi}}\times {s}_{{bi}}\right\rangle \left\langle {s}_{{di}}\times {s}_{{di}}\right\rangle }{\left\langle {s}_{b0}\times {s}_{b0}\right\rangle \left\langle {s}_{d0}\times {s}_{d0}\right\rangle }\right]}^{1/2}\\ \ \ +\,\mathrm{Cov}\left({n}_{a}\times {n}_{b},{n}_{c}\times {n}_{d}\right).\end{array}\end{eqnarray} \tag{ B3 }$$

The factors of 1/2 in the second and third terms are necessary to account for the factor of 2 difference between the variance of an auto-spectrum and the variance of the cross-spectrum between two uncorrelated fields.

An important feature of this rescaling method is that it considers signal and noise separately, rather than rescaling a combined signal-plus-noise covariance matrix. This is important because, for the highly nonuniform hit patterns (see Figure 7) that result from a large field-of-view telescope making deep maps, the spatial distribution of signal and noise in the map are very different. This leads to a significant difference in the number of signal versus noise degrees of freedom, and therefore different amounts of bandpower variance for the same overall power. This rescaling method does not attempt to account for different degrees of freedom between the different signal types, but that is a much smaller effect.