title: CMB-S4: Forecasting Constraints on Primordial Gravitational Waves
authors: CMB-S4 Collaboration : Kevork Abazajian Graeme E. Addison Peter Adshead Zeeshan Ahmed Daniel Akerib Aamir Ali Steven W. Allen David Alonso Marcelo Alvarez Mustafa A. Amin Adam Anderson Kam S. Arnold Peter Ashton Carlo Baccigalupi Debbie Bard Denis Barkats Darcy Barron Peter S. Barry James G. Bartlett Ritoban Basu Thakur Nicholas Battaglia Rachel Bean Chris Bebek Amy N. Bender Bradford A. Benson Federico Bianchini Colin A. Bischoff Lindsey Bleem James J. Bock Sebastian Bocquet Kimberly K. Boddy J. Richard Bond Julian Borrill François R. Bouchet Thejs Brinckmann Michael L. Brown Sean Bryan Victor Buza Karen Byrum Carlos Hervias Caimapo Erminia Calabrese Victoria Calafut Robert Caldwell John E. Carlstrom Julien Carron Thomas Cecil Anthony Challinor Clarence L. Chang Yuji Chinone Hsiao-Mei Sherry Cho Asantha Cooray Will Coulton Thomas M. Crawford Abigail Crites Ari Cukierman Francis-Yan Cyr-Racine Tijmen de Haan Jacques Delabrouille Mark Devlin Eleonora Di Valentino Marion Dierickx Matt Dobbs Shannon Duff Jo Dunkley Cora Dvorkin Joseph Eimer Tucker Elleflot Josquin Errard Thomas Essinger-Hileman Giulio Fabbian Chang Feng Simone Ferraro Jeffrey P. Filippini Raphael Flauger Brenna Flaugher Aurelien A. Fraisse Andrei Frolov Nicholas Galitzki Patricio A. Gallardo Silvia Galli Ken Ganga Martina Gerbino Vera Gluscevic Neil Goeckner-Wald Daniel Green Daniel Grin Evan Grohs Riccardo Gualtieri Jon E. Gudmundsson Ian Gullett Nikhel Gupta Salman Habib Mark Halpern Nils W. Halverson Shaul Hanany Kathleen Harrington Masaya Hasegawa Matthew Hasselfield Masashi Hazumi Katrin Heitmann Shawn Henderson Brandon Hensley Charles Hill J. Colin Hill Renée Hlozek Shuay-Pwu Patty Ho Thuong Hoang Gil Holder William Holzapfel John Hood Johannes Hubmayr Kevin M. Huffenberger Howard Hui Kent Irwin Oliver Jeong Bradley R. Johnson William C. Jones Jae Hwan Kang Kirit S. Karkare Nobuhiko Katayama Reijo Keskitalo Theodore Kisner Lloyd Knox Brian J. Koopman Arthur Kosowsky John Kovac Ely D. Kovetz Steve Kuhlmann Chao-lin Kuo Akito Kusaka Anne Lähteenmäki Charles R. Lawrence Adrian T. Lee Antony Lewis Dale Li Eric Linder Marilena Loverde Amy Lowitz Phil Lubin Mathew S. Madhavacheril Adam Mantz Gabriela Marques Frederick Matsuda Philip Mauskopf Heather McCarrick Jeffrey McMahon P. Daniel Meerburg Jean-Baptiste Melin Felipe Menanteau Joel Meyers Marius Millea Joseph Mohr Lorenzo Moncelsi Maria Monzani Tony Mroczkowski Suvodip Mukherjee Johanna Nagy Toshiya Namikawa Federico Nati Tyler Natoli Laura Newburgh Michael D. Niemack Haruki Nishino Brian Nord Valentine Novosad Roger O'Brient Stephen Padin Steven Palladino Bruce Partridge Don Petravick Elena Pierpaoli Levon Pogosian Karthik Prabhu Clement Pryke Giuseppe Puglisi Benjamin Racine Alexandra Rahlin Mayuri Sathyanarayana Rao Marco Raveri Christian L. Reichardt Mathieu Remazeilles Graca Rocha Natalie A. Roe Anirban Roy John E. Ruhl Maria Salatino Benjamin Saliwanchik Emmanuel Schaan Alessandro Schillaci Benjamin Schmitt Marcel M. Schmittfull Douglas Scott Neelima Sehgal Sarah Shandera Blake D. Sherwin Erik Shirokoff Sara M. Simon Anze Slosar David Spergel Tyler St. Germaine Suzanne T. Staggs Antony Stark Glenn D. Starkman Radek Stompor Chris Stoughton Aritoki Suzuki Osamu Tajima Grant P. Teply Keith Thompson Ben Thorne Peter Timbie Maurizio Tomasi Matthieu Tristram Gregory Tucker Caterina Umiltà Alexander van Engelen Eve M. Vavagiakis Joaquin D. Vieira Abigail G. Vieregg Kasey Wagoner Benjamin Wallisch Gensheng Wang Scott Watson Ben Westbrook Nathan Whitehorn Edward J. Wollack W. L. Kimmy Wu Zhilei Xu H. Y. Eric Yang Siavash Yasini Volodymyr G. Yefremenko Ki Won Yoon Edward Young Cyndia Yu Andrea Zonca
text: 
  Submitted to ApJ – draft version August 26, 2020
  
  
  
  
  
  
  
  CMB-S4—the next-generation ground-based cosmic microwave background (CMB) experiment—is set to
  significantly advance the sensitivity of CMB measurements and enhance 
  our understanding of the origin and evolution of the Universe, 
  from the highest energies at the dawn of time through the growth of structure to the present day. 
  Among the 
  science cases pursued with CMB-S4, 
  the quest for detecting primordial gravitational waves is a central driver of the experimental design. 
  This work details the development of a forecasting framework that includes
  a power-spectrum-based semi-analytic projection tool, 
  targeted explicitly towards optimizing constraints on the tensor-to-scalar ratio, r, 
  in the presence of Galactic foregrounds and gravitational lensing of the CMB. 
  This framework is unique in its direct use of information from the achieved performance 
  of current Stage 2–3 CMB experiments to robustly forecast the science reach 
  of upcoming CMB-polarization endeavors. 
  The methodology allows for rapid iteration over experimental configurations 
  and offers a flexible way to optimize the design of future experiments 
  given a desired scientific goal. 
  To form a closed-loop process, we couple this semi-analytic tool 
  with map-based validation studies, which allow for the injection of additional complexity and verification of our forecasts with several independent analysis methods. 
  We document multiple rounds of forecasts for CMB-S4 using this process and
  the resulting establishment of the current reference design of the primordial gravitational-wave component of the Stage-4 experiment, optimized to achieve our science goals of detecting primordial gravitational waves for r > 0.003 at greater than 5σ, or, in the absence of a detection, of reaching an upper limit of r < 0.001 at 95% CL.
  
  
  
  
  
  Victor Buza
  vbuza@kicp.uchicago.edu
  
  CMB-S4: Forecasting Constraints on Primordial Gravitational Waves
      
  
  
  CMBSFourAffiliationCount
  
  UCIrvine
  JohnsHopkinsUniversity
  UniversityofIllinoisatUrbana-Champaign
  SLAC
  UCBerkeley
  StanfordUniversity
  OxfordUniversity
  LawrenceBerkeleyNationalLaboratory
  RiceUniversity
  Fermilab
  UCSanDiego
  SISSA
  IFPU
  INFN
  HarvardUniversityCFA
  UniversityofNewMexico
  ArgonneNationalLaboratory
  AstroParticleandCosmologyLaboratory
  Caltech
  CornellUniversity
  KICP
  UniversityofChicago
  UniversityofMelbourne
  UniversityofCincinnati
  JPL
  LMUMunich
  UTAustin
  CITA
  InstitutdAstrophysiquedeParis
  StonyBrookUniversity
  UniversityofManchester
  ArizonaStateUniversity
  FloridaStateUniversity
  CardiffUniversity
  DartmouthCollege
  UniversityofGeneva
  UniversityofSussex
  UniversityofCambridge
  KavliInstituteforCosmologyCambridge
  IoA
  UniversityofTokyo
  UniversityofGroningen
  UniversityofToronto
  Dunlap
  KEK
  CEASaclay
  USTCH
  UniversityofPennsylvania
  McGillUniversity
  NIST
  PrincetonUniversity
  HarvardUniversity
  NASAGoddardSpaceFlightCenter
  SimonFraserUniversity
  UniversityofSouthernCalifornia
  HaverfordCollege
  StockholmUniversity
  CaseWesternReserveUniversity
  UniversityofBritishColumbia
  UniversityofColoradoBoulder
  UniversityofMinnesota
  UniversityofMichigan
  FlatironInstitute
  Columbia
  UniversityofVirginia
  KavliIPMU
  UCDavis
  YaleUniversity
  UniversityofPittsburgh
  BenGurionUniversity
  AaltoUniversity
  UCSantaBarbara
  PerimiterInstitue
  SouthernMethodistUniversity
  BerkeleyCenterforCosmologicalPhysics
  EuropeanSouthernObservatory
  GRAPPA
  WashingtonUniversitySt.Louis
  McDonnellCenterfortheSpaceSciencesWashingtonUniversity
  UniversityofMilano-Bicocca
  CentredePhysiquedesParticulesdeMarseille
  KIPAC
  InstituteforAdvancedStudy
  PennsylvaniaStateUniversity
  BrookhavenNationalLaboratory
  KyotoUniversity
  UniversityofWisconsinMadison
  UniversitedegliStudidiMilan
  ICJLab
  BrownUniversity
  SyracuseUniversity
  UCLosAngeles
  MSU
  The CMB-S4 Collaboration:
  Kevork Abazajian,<ref>
  Graeme E. Addison,<ref>
  Peter Adshead,<ref>
  Zeeshan Ahmed,<ref>
  Daniel Akerib,<ref>
  Aamir Ali,<ref>
  Steven W. Allen,<ref>
  David Alonso,<ref>
  Marcelo Alvarez,<ref>,<ref>
  Mustafa A. Amin,<ref>
  Adam Anderson,<ref>
  Kam S. Arnold,<ref>
  Peter Ashton,<ref>
  Carlo Baccigalupi,<ref>,<ref>,<ref>
  Debbie Bard,<ref>
  Denis Barkats,<ref>
  Darcy Barron,<ref>
  Peter S. Barry,<ref>
  James G. Bartlett,<ref>
  Ritoban Basu Thakur,<ref>
  Nicholas Battaglia,<ref>
  Rachel Bean,<ref>
  Chris Bebek,<ref>
  Amy N. Bender,<ref>,<ref>
  Bradford A. Benson,<ref>,<ref>,<ref>
  Federico  Bianchini,<ref>
  Colin A. Bischoff,<ref>
  Lindsey Bleem,<ref>,<ref>
  James J. Bock,<ref>,<ref>
  Sebastian Bocquet,<ref>
  Kimberly K. Boddy,<ref>,<ref>
  J. Richard Bond,<ref>
  Julian Borrill,<ref>,<ref>
  François R. Bouchet,<ref>
  Thejs Brinckmann,<ref>
  Michael L. Brown,<ref>
  Sean Bryan,<ref>
  Victor Buza,<ref>,<ref>
  Karen Byrum,<ref>
  Carlos Hervias  Caimapo,<ref>
  Erminia Calabrese,<ref>
  Victoria  Calafut,<ref>
  Robert Caldwell,<ref>
  John E. Carlstrom,<ref>,<ref>,<ref>
  Julien Carron,<ref>,<ref>
  Thomas Cecil,<ref>
  Anthony Challinor,<ref>,<ref>,<ref>
  Clarence L. Chang,<ref>,<ref>,<ref>
  Yuji Chinone,<ref>
  Hsiao-Mei Sherry Cho,<ref>
  Asantha Cooray,<ref>
  Will Coulton,<ref>
  Thomas M. Crawford,<ref>,<ref>
  Abigail Crites,<ref>,<ref>,<ref>
  Ari  Cukierman,<ref>,<ref>
  Francis-Yan Cyr-Racine,<ref>
  Tijmen de Haan,<ref>
  Jacques Delabrouille,<ref>,<ref>,<ref>
  Mark Devlin,<ref>
  Eleonora Di Valentino,<ref>
  Marion Dierickx,<ref>
  Matt Dobbs,<ref>
  Shannon Duff,<ref>
  Jo Dunkley,<ref>
  Cora Dvorkin,<ref>
  Joseph Eimer,<ref>
  Tucker Elleflot,<ref>
  Josquin Errard,<ref>
  Thomas Essinger-Hileman,<ref>
  Giulio Fabbian,<ref>
  Chang Feng,<ref>
  Simone Ferraro,<ref>
  Jeffrey P. Filippini,<ref>
  Raphael Flauger,<ref>
  Brenna Flaugher,<ref>
  Aurelien A. Fraisse,<ref>
  Andrei Frolov,<ref>
  Nicholas Galitzki,<ref>
  Patricio A. Gallardo,<ref>
  Silvia Galli,<ref>
  Ken Ganga,<ref>
  Martina Gerbino,<ref>
  Vera Gluscevic,<ref>,<ref>
  Neil Goeckner-Wald,<ref>
  Daniel Green,<ref>
  Daniel Grin,<ref>
  Evan Grohs,<ref>
  Riccardo Gualtieri,<ref>
  Jon E. Gudmundsson,<ref>
  Ian Gullett,<ref>
  Nikhel  Gupta,<ref>
  Salman Habib,<ref>
  Mark Halpern,<ref>
  Nils W. Halverson,<ref>
  Shaul Hanany,<ref>
  Kathleen Harrington,<ref>
  Masaya Hasegawa,<ref>
  Matthew Hasselfield,<ref>
  Masashi Hazumi,<ref>
  Katrin Heitmann,<ref>
  Shawn Henderson,<ref>
  Brandon Hensley,<ref>
  Charles Hill,<ref>
  J. Colin Hill,<ref>
  Renée Hložek,<ref>,<ref>
  Shuay-Pwu Patty Ho,<ref>
  Thuong Hoang,<ref>
  Gil Holder,<ref>
  William Holzapfel,<ref>
  John Hood,<ref>,<ref>
  Johannes Hubmayr,<ref>
  Kevin M. Huffenberger,<ref>
  Howard Hui,<ref>
  Kent Irwin,<ref>
  Oliver Jeong,<ref>
  Bradley R. Johnson,<ref>
  William C. Jones,<ref>
  Jae Hwan  Kang,<ref>
  Kirit S. Karkare,<ref>,<ref>
  Nobuhiko Katayama,<ref>
  Reijo Keskitalo,<ref>,<ref>
  Theodore Kisner,<ref>,<ref>
  Lloyd Knox,<ref>
  Brian J. Koopman,<ref>
  Arthur Kosowsky,<ref>
  John Kovac,<ref>,<ref>
  Ely D. Kovetz,<ref>
  Steve Kuhlmann,<ref>
  Chao-lin Kuo,<ref>
  Akito Kusaka,<ref>
  Anne Lähteenmäki,<ref>
  Charles R. Lawrence,<ref>
  Adrian T. Lee,<ref>,<ref>
  Antony Lewis,<ref>
  Dale Li,<ref>
  Eric Linder,<ref>
  Marilena Loverde,<ref>
  Amy Lowitz,<ref>,<ref>
  Phil Lubin,<ref>
  Mathew S. Madhavacheril,<ref>
  Adam Mantz,<ref>
  Gabriela Marques,<ref>
  Frederick Matsuda,<ref>
  Philip Mauskopf,<ref>
  Heather McCarrick,<ref>
  Jeffrey McMahon,<ref>,<ref>
  P. Daniel Meerburg,<ref>
  Jean-Baptiste Melin,<ref>
  Felipe Menanteau,<ref>
  Joel Meyers,<ref>
  Marius Millea,<ref>
  Joseph Mohr,<ref>
  Lorenzo Moncelsi,<ref>
  Maria Monzani,<ref>
  Tony Mroczkowski,<ref>
  Suvodip Mukherjee,<ref>,<ref>
  Johanna Nagy,<ref>,<ref>
  Toshiya  Namikawa,<ref>
  Federico Nati,<ref>
  Tyler Natoli,<ref>,<ref>
  Laura Newburgh,<ref>
  Michael D. Niemack,<ref>
  Haruki Nishino,<ref>
  Brian Nord,<ref>
  Valentine Novosad,<ref>
  Roger O'Brient,<ref>,<ref>
  Stephen Padin,<ref>
  Steven Palladino,<ref>
  Bruce Partridge,<ref>
  Don Petravick,<ref>
  Elena Pierpaoli,<ref>
  Levon Pogosian,<ref>
  Karthik Prabhu,<ref>
  Clement Pryke,<ref>
  Giuseppe Puglisi,<ref>,<ref>
  Benjamin Racine,<ref>,<ref>
  Alexandra Rahlin,<ref>
  Mayuri Sathyanarayana Rao,<ref>
  Marco Raveri,<ref>
  Christian L. Reichardt,<ref>
  Mathieu Remazeilles,<ref>
  Graca Rocha,<ref>
  Natalie A. Roe,<ref>
  Anirban Roy,<ref>
  John E. Ruhl,<ref>
  Maria Salatino,<ref>,<ref>
  Benjamin Saliwanchik,<ref>
  Emmanuel Schaan,<ref>
  Alessandro Schillaci,<ref>
  Benjamin Schmitt,<ref>
  Marcel M. Schmittfull,<ref>
  Douglas Scott,<ref>
  Neelima Sehgal,<ref>
  Sarah Shandera,<ref>
  Blake D. Sherwin,<ref>
  Erik Shirokoff,<ref>,<ref>
  Sara M. Simon,<ref>
  Anže Slosar,<ref>
  David Spergel,<ref>
  Tyler St. Germaine,<ref>
  Suzanne T. Staggs,<ref>
  Antony Stark,<ref>
  Glenn D. Starkman,<ref>
  Radek Stompor,<ref>
  Chris Stoughton,<ref>
  Aritoki Suzuki,<ref>
  Osamu Tajima,<ref>
  Grant P. Teply,<ref>
  Keith Thompson,<ref>
  Ben Thorne,<ref>
  Peter Timbie,<ref>
  Maurizio Tomasi,<ref>
  Matthieu Tristram,<ref>
  Gregory Tucker,<ref>
  Caterina Umiltà,<ref>
  Alexander van Engelen,<ref>
  Eve M. Vavagiakis,<ref>
  Joaquin D. Vieira,<ref>
  Abigail G. Vieregg,<ref>,<ref>
  Kasey Wagoner,<ref>
  Benjamin Wallisch,<ref>,<ref>
  Gensheng Wang,<ref>
  Scott Watson,<ref>
  Ben Westbrook,<ref>
  Nathan Whitehorn,<ref>,<ref>
  Edward J. Wollack,<ref>
  W. L. Kimmy Wu,<ref>
  Zhilei Xu,<ref>
  H. Y. Eric Yang,<ref>
  Siavash Yasini,<ref>
  Volodymyr G. Yefremenko,<ref>
  Ki Won Yoon,<ref>
  Edward Young,<ref>
  Cyndia Yu,<ref>
  Andrea Zonca<ref>
  
  
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  <ref>
  
  
  
      December 30, 2023
  ===============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================
  
  
  
  
  
  § INTRODUCTION
  
  
  
  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 <cit.>. 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, parametrized 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 <cit.>, 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 <cit.> 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 <cit.>. In addition to the production of PGWs <cit.>, 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  collaboration <cit.>. 
  
  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 <cit.>, BICEP/Keck <cit.>, CLASS <cit.>, POLARBEAR/Simons Array <cit.>, and SPT <cit.>, with Simons Observatory to follow soon <cit.>. Additionally, there are current and future balloon and satellite missions such as SPIDER <cit.> and LiteBIRD <cit.>, which we expect to complement ground-based measurements. The current best constraints are r(k_s=0.05 Mpc^-1, n_t=0)<0.06 at 95% confidence <cit.>, 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 <cit.>.
  To achieve these goals requires a significant increase in sensitivity, from 2–4×10^4 detectors in Stage 3 experiments to roughly 5×10^5 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 towards 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^5∼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 towards optimism, many of these forecasts have not in practice been achieved once the data have been taken and analyzed: 
  ABS  <cit.>; BICEP/Keck  <cit.>; EBEX  <cit.>; Planck  <cit.>; POLARBEAR  <cit.>; QUIET  <cit.>; SPIDER  <cit.>; SPTpol  <cit.>.
  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 multi-year 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 which 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 <cit.>. 
  
  
  
   §.§ 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
  non-idealities, 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 <ref>, to tie the semi-analytic 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.
  
  The main steps describing this process are as follows.
  
         
    * Develop a semi-analytic power-spectrum-level forecast, assuming noise performance that is scaled from analyses of real experiments. 
         
    * Use this forecasting tool to optimize the allocation of detector effort across observing frequencies, determining certain “checkpoints” in survey 
         definition space.
         
    * Use the checkpoint configurations to create standardized, version-numbered map-based data challenges (DCs) for validation.
         
    * Estimate science parameters from the DC maps with independent component-separation analysis methods.
         
    * 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.
         
    * 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.
  
  
  
  
  
  
   §.§ CMB-S4 r Forecasting Work-flow and Evolution
  
  Developing the forecasting machinery for CMB-S4, and increasing the
  robustness and realism of its results, has been an ongoing,
  ever-evolving, group effort, performed under the auspices of the
  CMB-S4 r Forecasting Working Group. 
  This work has undergone three major iterations, with results presented
  in the CMB-S4 Science Book <cit.>, the CMB-S4 Concept Definition Task force (CDT) Report <cit.>, and most recently in the CMB-S4 Reference
  Design Report <cit.>.
  
  Each stage in the evolution of this framework has been collectively
  considered, and carefully documented in the CMB-S4 Simulation and Forecasting Logbook.[<https://cmb-s4.org/wiki/index.php/Simulation_and_Forecasting_Logbook>]
  The resulting baseline survey definitions have been translated to
  publicly available, version-numbered map-based Data Challenges (DCs).[<https://cmb-s4.org/wiki/index.php/Data_Challenges>]
  
  Thus far, two independent groups have participated in testing the strengths and shortcomings of different analysis methods on these simulations, as well as improving the simulations themselves, as described in Sections <ref> and
  <ref>.
  
  Given the difficulty in describing a continuously-evolving process, for this paper we present several snapshots of our forecasting work and describe the relevant details. Sections <ref> to <ref> represent a full pass through the forecasting loop, as presented in Figure <ref>. In Sections <ref> and <ref> we present the full semi-analytic forecasting framework and optimization process, as used in the CMB-S4 CDT Report. In Section <ref> we discuss the creation of the DC4 simulation suite, corresponding to the baseline presented in that report. Sections <ref> and <ref> describe two independent analysis methods and the results obtained from applying these methods to DC4. In Section <ref> we describe several recent updates to the framework and the resulting findings presented in the CMB-S4 Reference Design Report. We conclude this paper with a discussion of the path forward. 
  
  
  
  
  § SEMI-ANALYTIC FORECASTING FRAMEWORK
  
  
  For the CMB-S4 Science Book <cit.>, we developed a semi-analytic forecasting framework specifically targeted towards 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 <cit.>. This automatically builds into the forecast all real-world inefficiencies including (but not limited to): imperfect detector yield; non-uniform detector performance; read-out noise; observing inefficiency; losses due to timestream filtering; beam smoothing; and non-uniform 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 <cit.>. We validate this approach using map-based simulations as part of the forecasting loop and present detailed results in Section <ref>. Our confidence in the projections is grounded in the connection to achieved performance and published results. 
  
  Figure <ref> presents a schematic representation of the semi-analytic forecasting framework, identifying the user inputs, code modules, and outputs. The subsections that follow describe this framework in detail.
  
  
  
  
  
  
   §.§ Fisher Formalism
  
  
  
  Given a likelihood function of the form
  
      L(θ;d)∝exp[-1/2(d-μ(θ))^TΣ(θ)^-1(d-μ(θ))]/√(det(Σ(θ))),
  
  where d are the data bandpowers, θ are the theory parameters, and μ(θ) and Σ(θ) are the bandpower expectation values and the bandpower covariance matrix given the parameters, we can calculate the expectation value of the log-likelihood curvature, evaluated at the position of the best fit model: 
  
      F_ij=-⟨∂^2logL(θ;d)/∂θ_i∂θ_j⟩.
  
  This quantity, called the Fisher information matrix, measures how steeply the likelihood falls as we move away from the best-fit model, and F^-1 can be thought of
  as the best possible covariance matrix for the measurement errors on the parameters θ_i. It can be shown that √((F^-1)_ii) is the minimum obtainable standard deviation on the desired parameters <cit.>.
  
  Inserting Equation <ref> into Equation <ref> yields
  
      F_ij=∂μ^T/∂θ_iΣ^-1∂μ/∂θ_j+1/2Tr(Σ^-1∂Σ/∂θ_iΣ^-1∂Σ/∂θ_j).
  
  We then calculate our parameter constraints as
  
      σ_i=√((F^-1)_ii).
  
  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 <ref> identically zero. Equation <ref> 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 <ref>).
  
  
  
   §.§ Forecasting Inputs
  
  In this section we briefly describe the set of inputs to the forecasting code, shown as green boxes in Figure <ref>. We discuss the key inputs in more detail in Sections <ref> and  <ref>.
  
  
  
  Achieved Performance: The code takes signal and noise simulations of the BICEP/Keck datasets <cit.> and derives the bandpower covariance matrix and the ensemble-averaged signal and noise bandpowers. These inputs contain information about the actual on-sky map noise achieved from multiple receivers, over multiple years, at 95, 150, and 220, including all real-world penalties. 
  A more detailed description of these simulations is available in Appendix H of <cit.>. Similar information from another experiment could easily be substituted. To project the performance of CMB-S4 channels, we assume that we can scale down the achieved noise based on increased detector count and integration time and that we can apply beam-size and noise-equivalent temperature (NET) rescalings to account for the differences in experimental design. 
  
  
  
  
  
  Scalable Instrument Specification: To specify the forecast instrument we start by selecting a set of observing frequency bands, with bandpass functions describing the response of each band. Then, for each observing band, we must provide the beam size, number of detectors, and ideal per-detector NET. The last two items can be used to make an idealized calculation of the instrument sensitivity in each band. We emphasize that we use these ideal performance numbers only for scaling between frequency bands, by comparing to similarly-calculated ideal sensitivities of BICEP/Keck. The scaling factors are ultimately applied to the achieved sensitivities, as described in Section <ref>, to obtain performance-based sensitivities in our desired bands. We note that for the frequency bands in which we do not currently have existing data, we extrapolate from the closest frequency for which we do.
  
  
  
  Fiducial Sky Model: A fiducial parametric model describing the foreground and CMB signal. Our standard model, which has 14 parameters, is discussed in Section <ref> and in Appendix A.
  
  
  
  
  Priors: If we have external prior knowledge for a given parameter θ_i, we can introduce this information by adding P_i=1/σ_i^2 to the diagonal of the Fisher matrix, where σ_i is the width of the prior. 
  
  
  
   §.§ CMB-S4 Scalable Instrument Specification
  
  
  
  
  To span the four atmospheric windows available to ground-based experiments (Figure <ref>) and have enough channels to mitigate against complex foregrounds, we assume eight channels at 30, 40, 85, 95, 145, 155, 220, and 270, 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 <cit.>, yielding the beams presented in Table <ref>. In addition, we also include a 20-GHz channel on a large-aperture telescope (LAT).[The inclusion of a 20-GHz channel is the result of insight gained from an earlier iteration through the forecasting loop (performed for the CMB-S4 Science Book), which demonstrated that for specific foreground models, sizeable biases were present due to synchrotron residuals. To mitigate against such biases the reference design was updated to include this additional channel.
  Placing this low-frequency band on a SAT would result in a very broad beam, which would dominate the noise at the relevant scales; to circumvent this, as mentioned above, we place it on a LAT. This means that while the scaling of the noise for this channel is still calculated from achieved performance, we use noise parameters (see Equation <ref>) that are more in tune with what has been achieved by LATs at the time of this forecasting <cit.>, i.e., the 1/f noise component is characterized by an ℓ_knee=200, keeping γ – the slope of this component – the same as for the small-aperture noise, and a beam of Θ_ν=11' FWHM (assuming a 6-m aperture).
  ]
  
  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[<https://cmb-s4.org/wiki/index.php/New_NET_Calculator_and_Validation>],[<https://github.com/dbarkats/NET_forecast_python>] at the South Pole and Chile, using the 10-year MERRA2 median atmospheric profiles <cit.>.
  We use the average of NET calculated for detectors at the two sites, which are 214, 177, 224, 270, 238, 309, 331, 747 and 1281K_CMB√(s) for our nine channels at 20–270, respectively.[Going beyond the number of bands for the current reference design increases manufacturing complexity and costs, and decreases the per-channel NET, thereby reducing the overall sensitivity in each channel; while exploring alternative options with 5 bands (20, 30, 95, 155, 270 GHz) and 7 bands (20, 30, 95, 155, 220, 270 GHz) has shown that choosing fewer bands leads to statistically significant biases on r after marginalization over foreground residuals (an increase in bias of up to Δ r=1.5×10^-4 when compared to the current count). Balancing these considerations, we have chosen the proposed configuration with nine frequency bands for this work. Going forward, as instrumentation choices are finalized, we anticipate a possible revision of this design.].
  These NETs are lower than similar BICEP/Keck channels (which are 287.6, 313.1 and 837.7K_CMB√(s) for 95, 150 and 220GHz) predominantly because they are calculated for a 100-mK thermal bath, as opposed to 250mK (which was also used in the CMB-S4 Science Book forecasts). 
  
  
  
  
  
  We also fold in information from two WMAP channels, 23 and 33, and seven Planck channels, 30, 44, 70, 100, 143, 217, and 353, 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.
  
  
  
    §.§.§ 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 <cit.>. Unlike Galactic foregrounds, the lensing signal is achromatic and cannot be mitigated with multi-frequency 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 <cit.>. 
  
  Reconstructing ϕ with high S/N requires high-sensitivity, high-angular resolution CMB polarization maps <cit.>. 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-arcminute 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 <ref>. 
  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  <ref> and Equation <ref>), but using SPTpol achieved performance <cit.>, i.e., without incurring penalties from non-idealities specific to low-resolution instruments and low-ℓ analysis, such as low-ℓ mode filtering and non-uniform coverage. Following the iterative formalism in <cit.>[Delensing estimators which are technically more optimal have been introduced in <cit.> and <cit.>, and we are currently developing tools to test their feasibility in further iterations of our forecasting.], using ℓ_min=300 and ℓ_max=4000 for the ϕ reconstruction and ℓ_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 <ref> and <ref>.
  
  
  
   §.§ Multicomponent Theory Model
  
  
  
  
  Our model includes a CMB component parametrized by r and the residual lensing amplitude, A_L, and components of polarized dust and synchrotron emission. We assume that the synchrotron scales as a simple power law in both frequency and ℓ. For the dust we assume a power-law scaling in ℓ and a modified blackbody spectral energy distribution (SED). We allow for spatially correlated synchrotron and dust, parameterized by a single correlation parameter ϵ <cit.>; the effective frequency scaling of this correlation depends on the relative strength of the two components. In addition, we also consider dust and synchrotron frequency decorrelation parameters, which allow their spatial pattern to change with frequency, suppressing the correlation of foreground signals between observing bands. A detailed description of the full parametric model is presented in Appendix <ref>. The current model is easily extendable to accommodate additional complexities that have not yet been captured; alternatively, other models could be substituted here as well.
  
  
  
  The model parameters are:
  
    
  
    * r, tensor-to-scalar ratio, at pivot scale k_s=0.05  Mpc^-1
  
       and spectral index of the tensor modes n_t=0;
  
    * A_L, residual lensing amplitude;
  
    * A_dust, dust amplitude, in μK^2_CMB, 
  
           at 353 GHz and ℓ=80;
  
    * β_d, dust spectral index;
  
    * T_d, dust temperature;
  
    * α_d, dust spatial spectral index;
  
    * Δ_d, dust frequency correlation, 
  
         between 217 and 353 GHz, at ℓ=80;
  
    * EE/BB_dust, power ratio for dust;
  
    * A_sync, synchrotron amplitude, in μK^2_CMB, 
  
           at 23 GHz and ℓ=80;
  
    * β_s, synchrotron spectral index;
  
    * α_s, synchrotron spatial spectral index;
  
    * Δ_s, synchrotron frequency correlation,
  
         between 23 and 33 GHz, at ℓ=80;
  
    * EE/BB_sync, power ratio for synchroton;
  
    * ϵ, dust/synchrotron spatial correlation.
  
  
  For a given auto- or cross-spectrum, we step through the model components, combine the appropriate amplitude functions for the two bands contributing to the spectrum, and apply the bandpower window functions to obtain the binned expectation values. Finally, we sum over model components to find the total expectation value for that spectrum. In addition, since a Fisher forecast requires knowledge of the response of the model expectation values with respect to the model parameters, we also output the derivatives of the model expectation values.
  
  The Fisher matrix that we consider is usually 10-dimensional. The parameters we constrain are r, A_dust, β_d, α_d, Δ_d, A_sync, β_s, α_s, Δ_s, ϵ. We fix
  T_d=19.6 because this parameter is mostly degenerate with A_dust for observations below 300, where the SED is in the Rayleigh-Jeans limit. The parameter A_L is assumed to be known, but its value is adjusted to represent varying levels of delensing, as discussed in Section <ref>. The
  EE/BB ratios are not relevant for calculations presented here because we are focusing on constraints from the BB spectrum only; however, these parameters are left in for possible future forecasting. 
  The fiducial model used for forecasting is centered at either r = 0 or 0.003, with A_dust=4.25 μK^2_CMB <cit.> and A_sync=3.8 μK^2_CMB <cit.>. In Sections <ref>–<ref> we assume no variation of these foreground amplitudes over the sky fractions relevant to this study, i.e., they are always pinned to the values listed above. We revisit this assumption in Section <ref>. The spatial and  frequency spectral indices are centered at the preferred Planck and WMAP values <cit.>: 
  β_d=1.59 (with Gaussian prior of width 0.11); β_s=-3.10 (with a Gaussian prior of width 0.30); α_d=-0.42; α_s=-0.6; and the dust/synchrotron correlation is centered at ϵ=0. Following <cit.>, the central dust correlation value is taken to be Δ_d=0.97 (3% decorrelation) and 
  the synchrotron correlation value is assumed to be Δ_s=1 (no fiducial synchrotron decorrelation). Unless otherwise stated, the parameters have flat unbounded priors.  
  
  
  
   §.§ Bandpower Covariance Matrix Rescaling
  
  
  
  
  
  
    §.§.§ 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 Cov(𝒩_i×𝒩_j, 𝒩_k×𝒩_l);
    
    * sn1 = signal×noise terms Cov(S_i×𝒩_j,S_k×𝒩_l);
    
    * sn2 = signal×noise terms Cov(S_i×𝒩_j, 𝒩_k× S_l);
    
    * sn3 = signal×noise terms Cov(𝒩_i× S_j, S_k×𝒩_l);
    
    * sn4 = signal×noise terms Cov(𝒩_i× S_j, 𝒩_k× S_l).
  
  Here S are signal simulations, 𝒩 are 
  noise simulations, and the indices i, j, k, 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 <ref>. 
  
  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.
  
  
  
  
    §.§.§ 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 dataset 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 dataset and scale down the noise in the BPCM by the desired amount. In particular, each BPCM component is scaled independently by √(N_ℓ,S4/N_ℓ,BK) for each factor of 𝒩 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_ℓ^2. The survey weight is defined as w=2f_sky/σ_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
  
      N_ℓ,S4=N_ℓ,BKw_BK,achieved/w_S4,projectedB^2_ℓ,S4/B^2_ℓ,BK,
  
  where B^2_ℓ,ν=exp-ℓ(ℓ+1)Θ^2_ν/8log(2), Θ_ν is the full width at half maximum (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 
  
      w_S4,projected/w_BK,achieved=n_S4^det-yr/n_BK^det-yrNET^2_BK,ideal/NET^2_S4,ideal,
  
  where n_S4^det-yr is the number of detector-years assumed for CMB-S4 at any particular frequency, and n_BK^det-yr
  is the number of detector-years in the BK15 dataset, with the instruments in their final state <cit.>.
  
  
  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 are scaled from BICEP/Keck 150-GHz data; and 220 and 270 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^BK_sky=Ω_pix/4π∑_im_i≃ 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 β=f^S4_sky
  /f^BK_sky, which boosts the (signal × noise) and (noise × noise) terms of the covariance matrix by β and β^2, 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 non-uniform coverage have on the covariance structure. We revisit the way this scaling is performed in Section <ref>, 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.
  
  
  
  § OPTIMIZED FORECASTING FOR R
  
  
  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 (similar to a BICEP/Keck 150-GHz receiver observing for 1 year). For other channels, the number of detector-years per unit of effort is calculated as n_det,150×(ν/150)^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). 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^6 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 <ref> 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^-4, for f_sky=3% (see discussion on sky fraction below), requires 1.2× 10^6 150-GHz equivalent detector-years (or 1.8 × 10^6 when including marginalization over the decorrelation parameters). 
  
  Upon obtaining the optimized detector count distribution, we obtain the input noise spectra according to Equation <ref>. To use these spectra to create noise simulations (discussed in Section <ref>), it is useful to distill them to a few input parameters. To that end we fit them to the formula
  
      N_ℓ,fit = ℓ(ℓ+1)/2πΩ_pix/B_ℓ^2(1+(ℓ/ℓ_knee)^γ)σ_map^2
  
  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 <ref> and the configuration that achieves the science goal in 4 years (or 6 years when marginalizing over decorrelation) is summarized in Table <ref>.
  
  
  
  
  
  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 <ref> 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. 
  
  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 <ref>, 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 <cit.>.
  
  Second, the current optimization assumes uniform foreground behavior across the sky (with amplituide 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 <cit.>. 
  
  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 <cit.> and CMB-S4 CDT Report <cit.> versions of the forecasts. As mentioned in Section <ref>, we revisit the issue of sky coverage in Section <ref> with updated assumptions about our survey strategy and how we perform the f_sky scaling from the achieved products. 
  
  
  
    
  
  
  
  
  
  § MAP-BASED SIMULATIONS
  
  
  
  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 semi-analytic 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 <ref>.
  
  
  
   §.§ Map Noise Simulations
  
  
  We use Equation <ref> to obtain the desired noise prescription, for a fixed total effort, and then generate 499 Gaussian noise realizations at each band. As in the Science Book, we have mostly used a circular sky area of 3%. Small-aperture cameras have a very wide instantaneous field of view and hence the observed sky region necessarily has a large edge taper. For the nominal 3% sky coverage simulations, we assumed a circular sky patch centred at RA=0, Dec=-45 (slightly below the BICEP/Keck patch) with full coverage out to radius θ<12 and “relative hits” tapering to zero with a cosine-squared shape for 12<θ<27. The noise realizations are divided by the square root of this coverage pattern such that the noise “blows up around the edge” as it does in real maps. We also assume an ℓ_min=30 cut-off below which we do not recover any information.
  
  
  
  For the CMB-S4 Reference Design Report (and Section <ref>), we include an explicit scan strategy on the sky and produce more realistic sky coverage patterns, but for the moment we regard the above as a reasonable compromise between idealism and reality.
  
  
  
   §.§ Foreground Models
  
  
  
  To make simulated sky maps, we add realizations of lensed CMB (both with and without a PGW component) to models of the Galactic foregrounds.  So far, we have run simulations with seven foreground models, which we now describe.
  
  
  
  
  
    * Simple Gaussian realizations of synchrotron and dust with power-law angular power spectra at amplitudes set to match the observations in the BICEP/Keck field, and simple uniform SEDs (power law for synchrotron, modified blackbody for dust).
  
  
    * The PySM[<https://github.com/bthorne93/PySM_public>] model a1d1f1s1, where the letters refer to anomalous microwave emission (AME), dust, free-free, and synchrotron, respectively, and the numbers are the base models described in <cit.>. Free-free and AME are assumed to be unpolarized in this model
  and thus do not affect the analysis in this paper.
  
  
    * The PySM model a2d4f1s3, where the models have been updated to variants that are also described in <cit.>.  Note that these include 2% polarized AME, a curvature of the synchrotron SED, and a two-temperature model for dust.
  
  
    * The PySM model a2d7f1s3, where the dust model has been updated to a more sophisticated physical characterization of dust grains as described in <cit.>. This model is interesting in that it does not necessarily conform to the modified blackbody SED.
  
  
    * The dust in Model 3 is replaced by a model of polarized dust emission that incorporates Hi column density maps as tracers of the dust intensity structures, and a phenomenological description of the Galactic magnetic field as described in <cit.>.  The model is expanded beyond what is described in that paper to produce a modest amount of decorrelation of the dust emission pattern as a function of frequency motivated by the analysis of Planck data in <cit.>.
  
  
    * A toy model where the strong dust decorrelation suggested in Figure 3 of <cit.> is taken at face value (Δ^217×353=0.85, at ℓ=80) and scaled to other frequencies using the functional form given in appendix B of <cit.>, with a linear scaling in ℓ.
  While such a model is not ruled out by current data, it appears to be very hard to produce such strong decorrelation in physics-based models.  We also note that <cit.> and <cit.> have re-analyzed the same Planck data and, while they find that the high level of decorrelation in this model is still consistent with the data, their best fit to that same data has no decorrelation.
  
  
    * A model based on MHD simulations <cit.> of the Galactic magnetic field, which naturally produces non-Gaussian correlated dust and synchrotron emission.
  
  
  
  Models 1 to 4 use the large-scale modes of the real sky as measured above the noise in the Planck data.  This means that these models are intrinsically “single-realization,” and this must be borne in mind when interpreting the results. Models 4 and 6 are not based on Planck data, but still contain a fixed signal realization. Models 0 and 5 have different seeds for each signal map and include the (Gaussian) sample variance. The PySM models fill in the small-scale structure with power-law Gaussian extrapolations, while Models 4 and 6 naturally produce non-Gaussian small-scale structure. However, all of these models are consistent with current data, and the more complex models are not necessarily more accurate reflections of reality.
  
  
  
   §.§ 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. <cit.>,  <cit.>, <cit.>, <cit.>, <cit.>.
  There have also been several publications that examine the effect of specific classes of instrumental systematics on a generic polarization experiment, e.g. <cit.>, <cit.>, <cit.>, <cit.>, <cit.>.
  
  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 re-analyzing 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 semi-analytic 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.
  
  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.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Delensing
  
  
  
  We have generated high-resolution simulated maps on which we intend to run explicit lensing reconstruction and then include that information in the analysis. While we are currently working on this analysis, this process has not yet converged, and so for the present, we approximate delensing by scaling down the ΛCDM lensing signal by the appropriate factor, as described in Section  <ref>.
  
  
  
  
  § SIMULATION ANALYSIS METHODS
  
  
  
  To make simulated maps, the noise realizations described in Section <ref> are added to the sky models described in Section <ref>.
  For each realization, one then has a stack of multi-frequency I,Q,U maps containing non-uniform noise, foregrounds, and signal, and the challenge is to re-analyze 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 <cit.>, 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  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 multi-component fit to the ensemble of auto- and cross-spectra, as used for the BICEP/Keck analysis to date <cit.>.  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 <ref>.
  
  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.
  
  
  
  
  § MAP-BASED RESULTS
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Table <ref> summarizes the results of the analysis for simulations of the optimized configuration obtained in Section <ref>
  (1.2×10^6 150-GHz-equivalent detector-years) and residual lensing power
  A_L=0.1. The lensing residual is expected for iterative EB delensing according to <cit.> 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^-4, exactly as
  expected from the semi-analytic formalism.  As we progress to the more complex
  foreground models, σ(r) is generally in the range 5–8×10^-4. 
  
  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 approximatively 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 <ref> 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 years of operation is r=0.004.[The reason this is larger than 5 times the quoted σ(r) is sample variance.] For a tensor-to-scalar ratio of r=0.003, the median detection significance after 4 years is expected to be 4σ.  If a detection were to be emerging at this point, extending the run time to 8 years 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 years 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 <ref> 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^-4. The dominant contribution to the bias comes from synchrotron residuals, which motivated placing one lower-frequency channel on the LAT (reducing the magnitude of the bias to 1.8 × 10^-4). 
  
  
  
  Table <ref> 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^-4 in semi-analytic 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^-4 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 <cit.>.  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 parametrized models of the frequency spectra of different foregrounds, shows biases on r at the 1 × 10^-4 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.
  
  
  
  
  
  
  § MODIFICATIONS LEADING TO THE REFERENCE DESIGN
  
  
  
  In this section, we describe updates to the framework and the reference design, which 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 <cit.>. 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 <ref>. 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 <cit.> with 18 SAT tubes, observing for 7 years. We use this configuration to scale the BICEP/Keck noise bandpower covariance matrix in the same way as described in Section <ref>. A further re-scaling is applied to account for sky coverage, as explained below.
  
  
  
   §.§ Sky Coverage Effects
  
  The semi-analytic calculations of Section <ref>
  assumed a simplified re-scaling for sky area, while the map based simulations of Section <ref> assumed an idealized circular sky patch, which is not actually achievable with a practical instrument from a site at any latitude.
  Figure <ref> compares our prior assumptions to more realistic hit patterns.
  
  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 <ref> 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 Right Ascension. 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 <ref> have the effective sky fractions reported in Table <ref>, assuming inverse noise-variance weighting.
  
  
  
  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
  
      f_ eff^ noise =   Ω_ pix/4π∑_i  w_i^2 h_i^-1∑_i  h_i/∑_i w_i^2,
  
  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
  
  
      f_ sky^ noise =   Ω_ pix/4π(∑_i  w_i^2 h_i^-1)^2/∑_i w_i^4 h_i^-2,
  
  
      f_ sky^ signal =   Ω_ pix/4π(∑_i w_i^2)^2/∑_i w_i^4,
  
  
      f_ sky^ cross =  Ω_ pix/4π∑_i w_i^2 ∑_i  w_i^2 h_i^-1/∑_i w_i^4 h_i^-1.
  
  
  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.
  
  Lastly, 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  353 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.
  
  
  
   §.§ Delensing Revised
  
  The CMB-S4 PGW science goal can only be achieved if the majority of the lensing B modes can be removed. The optimization in Section <ref> assumed a single frequency channel assigned to the higher resolution delensing
  observations. The strength of polarized foregrounds at small angular scales is currently poorly constrained by data; the reference design, therefore, includes some additional coverage at higher and lower frequencies.[Though studies related to the role polarized small-scale galactic foreground residuals play in the delensing process have been performed in <cit.> and <cit.>, a complete understanding of their impact is still missing.]
  
  To forecast the delensing performance, we proceed in two steps. For a given LAT configuration and sky coverage, we derive the noise levels for an ILC that minimizes the variance of components with a frequency dependence that differs from that of a blackbody <cit.>. In this step, we assume that polarized foreground emission is dominated by Galactic synchrotron and thermal dust emission. Using the ILC noise power spectrum, we then forecast the performance expected for iterative EB delensing <cit.>.
  
  The LAT reference design, established independently from this work, and described in <cit.>, includes 2 LATs in Chile for the wide field survey targeting science goals other than PGWs, and 1 LAT at the South Pole for the PGW delensing survey. For the Chile LATs, and a wide-area survey covering 70% of the sky, the two-step procedure predicts that 73% of the lensing power can be removed in the “Chile shallow” region after 7 years of observation. Similarly, for the single LAT at the South Pole dedicated to delensing of the approximately 3% “Chile deep” and “Pole deep” regions, we expect to be able to remove close to 90% of the lensing power after 7 years of observation.
  
  
  
  
  The numbers given above assume identical hits maps for the South Pole LATs and SATs, and assume inverse noise-variance weighting rather than a weighting scheme that accounts for both signal and noise. For the lensing residuals achieved by the reference design, inverse noise-variance weighting for the SAT maps is suboptimal because the signal (e.g., for r=0 lensing residual after foreground removal) is measured with signal-to-noise ratio above unity over a region that extends beyond the region that dominates the analysis in inverse noise-variance weighting. As a consequence, in all the forecasts presented below we employ weights that account for both signal and noise, which significantly increases the number of modes that contribute to the analysis, especially for small, deep patches. Even though, for the same survey, this leads to slightly higher noise and lensing residuals, the increase in the number of modes overall reduces the lensing sample variance contribution to σ(r). In practice we determine the optimal weights iteratively, accounting for the effect on the SAT and LAT analyses and find rapid convergence.
  
  
  
   §.§ 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 <ref> 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 sub-region using the appropriate delensing level for that sub-region. 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.
  
  We also explore the possibility of unmodeled foreground residuals contributing residual power to the cleaned maps, assuming a foreground cleaning to 1% at 95GHz. 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 95GHz 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 <ref>, 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 <ref>). 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.
  
  In Tables <ref>–<ref> we present a set of σ(r) results for seven years of observations and r=0, and in Tables <ref>–<ref> 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.
  
  
  
  
  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 <ref>), there is still a difference
  between σ(r) obtained for equivalent numbers of optics tubes in Chile versus 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 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 <ref>.
  
  
  
  § CONCLUSIONS
  
  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 semi-analytic calculations with map-based simulations 
  of increasing complexity, thereby creating a closed forecasting loop (presented in Figure <ref>). 
  
  
  
  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 semi-analytic calculations described in Sections <ref> and <ref> indicate that for a 3% sky fraction, 1.2× 10^6 150-GHz-equivalent detector-years (or 1.8 × 10^6 when including marginalization over the decorrelation parameters) are required to reach the science requirement of σ(r)=5×10^-4, with 30% of this effort assigned to the delensing survey. The resulting optimal distribution of detectors across frequency bands is summarized in Table <ref>.
  
  The map-based simulations described in Sections <ref>–<ref> confirm the σ(r) results from the semi-analytic 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 semi-analytic calculations onto realizable instruments (see Section <ref>) results in the reference design described in <cit.> for a seven-year survey period, and in the constraints presented in Tables <ref>–<ref>.
  
  In Section <ref> we have also extended the semi-analytic 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 parametrization; 
  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.
  
  
  
  § ACKNOWLEDGEMENTS
  
  
  
  
  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 US, work on CMB-S4 by individual investigators has been supported by the National Science Foundation (awards 1248097, 1255358, 1815887, 1835865, 1852617), the Department of Energy (awards DE-SC0009919, DE-SC0009946, DE-SC0010129), 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, RH is supported by the Discovery Grants program from NSERC, and acknowledges funding from CIFAR, the Sloan Foundation and the Dunlap family.
  In Italy, CB acknowledges support under the ASI COSMOS and INFN INDARK programs.
  In the Netherlands, DM acknowledges NWO VIDI award number 639.042.730.
  In Switzerland, JC is supported by a SNSF Eccellenza Professorial Fellowship (No. 186879).
  In the United Kingdom, AL, GF, JC are supported by the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement No. [616170].  AL also acknowledges STFC award ST/P000525/1.
  SM 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 programme 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.
  
  
  
  
  
  
  
  
  § MULTICOMPONENT THEORY MODEL
  
  
  
  
  
  The forecasting framework uses a parametric model to describe the bandpower expectation values as a combination of cosmological and foreground signals. This parametrization follows the one presented in Appendix G of <cit.>.
  
  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 ν_1 and ν_2 can be written in the following form, after computing the expected spectrum as a function of ℓ, and applying bandpower window functions are applied to calculate binned bandpower expectation values:
  
      𝒟_ℓ,BB^ν_1 ×ν_2 = 
        A_dust f_d^ν_1 f_d^ν_2( ℓ/80)^ + 
        Δ'_s f_s^ν_1 f_s^ν_2( ℓ/80)^ 
        + ϵ√(A_dust) (f_d^ν_1 f_s^ν_2 + f_s^ν_1 f_d^ν_2) ( ℓ/80)^( + ) / 2.
  
  
  
  
  
  
  
  
  
  
  Parameters A_dust and  specify the dust and synchrotron power in units of μK_^2 at angular scale ℓ = 80.
  These are defined at pivot frequencies of 353 for dust and 23 for synchrotron.
  The dust and synchrotron components scale as power laws in ℓ with slopes  and , respectively (note that slope parameters are defined for 𝒟_ℓ≡).
  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 re-analysis steps. 
  
  Coefficients f_d^ν and f_s^ν, given by equations <ref> and <ref>, describe the scaling of dust and synchrotron amplitude from the pivot frequencies to the actual bandpasses of the maps at frequencies ν_1 and ν_2. 
  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  <cit.>.  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  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 <cit.>, 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
  
      f_MBBS^ν = ∫dν R(ν) ν^3+β_d(exphν/kT_d - 1)^-1/ν_pivot^3+β_d(exphν_pivot/kT_d - 1)^-1 .
  
  For a synchrotron power-law scaling (PLS) with spectral index β_s, we calculate the coefficient in a similar way,
  
      f_PLS^ν = ∫dν R(ν) ν^2+β_s/ν_pivot^2+β_s .
  
  The conversion between μK_ units at the pivot frequency and μK_ units at the target map bandpass is given by
  
      f_T_CMB^ν = ν_pivot^4 exphν_pivot/kT_CMB(exphν_pivot/kT_CMB-1)^-2/∫dν R(ν) ν^4 exphν/kT_CMB(exphν/kT_CMB-1)^-2 .
  
  Combining these factors, we obtain the scalings used in equation <ref>, 
  
      f_d^ν   = f_MBBS^ν× f_T_CMB^ν
      
        f_s^ν   = f_PLS^ν× f_T_CMB^ν
  
  
  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 <cit.>. 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
  
      =
      𝒟_80(217× 353)/√(𝒟_80(217× 217)𝒟_80(353× 353)) ,
  
  where 𝒟_80 is the dust power at ℓ=80.
  Here <1 corresponds to decorrelation. We scale to other frequency combinations using the factor suggested by <cit.>,
  
      f(ν_1,ν_2) = (log (ν_1/ν_2))^2/(log(217/353))^2.
  
  Similarly, based on suggestions from <cit.>, we consider three possible scalings with ℓ:
  
      g(ℓ) = {[              1,      flat case;;       (ℓ / 80),    linear case;;     (ℓ / 80)^2, quadratic case. ].
  
  
  The scalings above can produce extreme (and non-physical)
  behavior at high ℓ or for frequencies that are widely separated. We therefore re-map the scaled correlation coefficient using the following function:
  
      (ν_1,ν_2,ℓ) = exp[ log()   f(ν_1,ν_2)   g(ℓ) ] .
  
  With this re-mapping,  remains in the
  range 0 (no correlation) to 1 (perfect correlation) for all values of f and g.
  This combination of frequency scaling and non-linear re-mapping has been shown to correspond to a Gaussian spatial variation in the foreground spectral index parameter <cit.>.
  
  In a similar vein, we define the parameter Δ'_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
  
      D_ℓ,BB^ν_1×ν_2=r/0.1 D_ℓ,BB^tensor+A_L D_ℓ,BB^lensing,
  
  where D_ℓ,BB^tensor is the BB template for a tensor signal with r=0.1 and D_ℓ,BB^lensing is the expected lensing BB spectrum for ΛCDM. These are obtained using the [<https://camb.info/>] package <cit.>.
  
  
  
  
  
  
  
  
  § BPCM CONSTRUCTION AND RESCALING
  
  
  
  The bandpower covariance matrix (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 physically unrelated signals are 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 <cit.>, 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:  ⟨ s_xi× s_yj⟩ = 0 for i  j (independent signal fields)
  and ⟨ s_xi× n_y ⟩ = 0 (signal is independent of noise). Then, a generic bandpower covariance term can be written as
  
      Cov( m_a × m_b, m_c × m_d )    = ⟨ ( m_a × m_b ) ( m_c × m_d ) ⟩ - ⟨ m_a × m_b ⟩⟨ m_c × m_d ⟩
         = ∑_i Cov( s_ai× s_bi, s_ci× s_di) 
            + ∑_i ∑_j  iCov( s_ai× s_bj, s_ci× s_dj) + ∑_i ∑_j  iCov( s_ai× s_bj, s_cj× s_dj) 
            + ∑_i Cov( s_ai× n_b, s_ci× n_d ) + ∑_i Cov( s_ai× n_b, n_c × s_di) 
            + ∑_i Cov( n_a × s_bi, s_ci× n_d ) + ∑_i Cov( n_a × s_bi, n_c × s_di) 
            + Cov( n_a × n_b, n_c × n_d ) ,
  
  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., ⟨ n_x × n_y ⟩ = 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
  
      Cov( m_a × m_b, m_c × m_d )    = Cov( s_a0× s_b0, s_c0× s_d0) ∑_i [ ⟨ s_ai× s_ai⟩⟨ s_bi× s_bi⟩⟨ s_ci× s_ci⟩⟨ s_di× s_di⟩/⟨ s_a0× s_a0⟩⟨ s_b0× s_b0⟩⟨ s_c0× s_c0⟩⟨ s_d0× s_d0⟩]^1/2
         + 1/2Cov( s_a0× s_b0, s_c0× s_d0) ∑_i ∑_j  i[ ⟨ s_ai× s_ai⟩⟨ s_bj× s_bj⟩⟨ s_ci× s_ci⟩⟨ s_dj× s_dj⟩/⟨ s_a0× s_a0⟩⟨ s_b0× s_b0⟩⟨ s_c0× s_c0⟩⟨ s_d0× s_d0⟩]^1/2
         + 1/2Cov( s_a0× s_b0, s_c0× s_d0) ∑_i ∑_j  i[ ⟨ s_ai× s_ai⟩⟨ s_bj× s_bj⟩⟨ s_cj× s_cj⟩⟨ s_di× s_di⟩/⟨ s_a0× s_a0⟩⟨ s_b0× s_b0⟩⟨ s_c0× s_c0⟩⟨ s_d0× s_d0⟩]^1/2
         + Cov( s_a0× n_b, s_c0× n_d ) ∑_i [ ⟨ s_ai× s_ai⟩⟨ s_ci× s_ci⟩/⟨ s_a0× s_a0⟩⟨ s_c0× s_c0⟩]^1/2
         + Cov( s_a0× n_b, n_c × s_d0) ∑_i [ ⟨ s_ai× s_ai⟩⟨ s_di× s_di⟩/⟨ s_a0× s_a0⟩⟨ s_d0× s_d0⟩]^1/2
         + Cov( n_a × s_b0, s_c0× n_d ) ∑_i [ ⟨ s_bi× s_bi⟩⟨ s_ci× s_ci⟩/⟨ s_b0× s_b0⟩⟨ s_c0× s_c0⟩]^1/2
         + Cov( n_a × s_b0, n_c × s_d0) ∑_i [ ⟨ s_bi× s_bi⟩⟨ s_di× s_di⟩/⟨ s_b0× s_b0⟩⟨ s_d0× s_d0⟩]^1/2
         + Cov( n_a × n_b, n_c × n_d ) .
  
  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 non-uniform hit patterns (see Figure <ref>) 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.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  