# Cosmological Parameter Forecasts for a CMB-HD Survey

Amanda MacInnis  Physics and Astronomy Department, Stony Brook University,
Stony Brook, NY 11794  Neelima Sehgal  Physics and Astronomy Department, Stony
Brook University, Stony Brook, NY 11794  Miriam Rothermel  Physics and
Astronomy Department, Stony Brook University, Stony Brook, NY 11794

###### Abstract

We present forecasts on cosmological parameters for a CMB-HD survey. For a
$\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$ model, we find
$\sigma(n_{\mathrm{s}})=0.0013$ and $\sigma(N_{\mathrm{eff}})=0.014$ using CMB
and CMB lensing multipoles in the range of $\ell\in[30,20000]$, after adding
anticipated residual foregrounds, delensing the acoustic peaks, and adding
DESI BAO data. This is about a factor of two improvement in ability to probe
inflation via $n_{\mathrm{s}}$ compared to precursor CMB surveys. The
$N_{\mathrm{eff}}$ constraint can rule out light thermal particles back to the
end of inflation with 95% CL; for example, it can rule out the QCD axion in a
model-independent way assuming the Universe’s reheating temperature was high
enough that the axion thermalized. We find that delensing the acoustic peaks
and adding DESI BAO tightens parameter constraints. We also find that baryonic
effects can bias parameters if not marginalized over, and that uncertainties
in baryonic effects can increase parameter error bars; however, the latter can
be mitigated by including information about baryonic effects from kinetic and
thermal Sunyaev-Zel’dovich measurements by CMB-HD. The CMB-HD likelihood and
Fisher estimation codes used here are publicly available; the likelihood is
integrated with Cobaya to facilitate parameter forecasting.

##  I Introduction

Measurements of the Cosmic Microwave Background (CMB) have provided precise
information on the cosmological parameters of our Universe [1, 2, 3]. Near
future CMB experiments such as the Simons Observatory (SO) [4] and CMB-S4 [5]
will improve measurements of the CMB temperature and polarization power
spectra below $\ell\sim 5000$ and the CMB lensing power spectrum considerably,
yielding improved constraints on $\Lambda$CDM parameters, the sum of the
neutrino masses ($\sum m_{\nu}$), and the effective number of light
relativistic species ($N_{\mathrm{eff}}$), in particular [4, 5]. Beyond that,
an even lower noise and higher resolution CMB experiment can improve parameter
constraints further by measuring smaller angular scales in both temperature
and polarization, measuring the lensing power spectrum with more precision and
over smaller scales, and enabling more aggressive foreground removal down to
lower flux limits. CMB-HD is a proposed concept for a Stage-V CMB facility [6,
7] that would have three times lower instrument noise and about six times
higher resolution than CMB-S4. In this work, we forecast the cosmological
parameter constraints that can be achieved by such a facility, folding in
anticipated residual foregrounds, delensing, and baryonic effects.

Delensing of the CMB acoustic peaks has been demonstrated in a number of
recent CMB data analyses [8, 9, 10]. It can tighten cosmological parameter
constraints by sharpening the acoustic peaks, removing uncertainty in the
lensing realization, and reducing lensing-induced power spectrum covariances
[11, 12]. Baryonic effects such as AGN feedback can move around the matter
distribution on moderately small scales, which impacts both the lensed CMB
power spectra and the CMB lensing power spectrum; for next-generation CMB
experiments, neglecting such effects can bias parameters [13, 14]. The
residual foregrounds we include in this work, in both CMB and CMB lensing
power spectra, are those that preliminary analyses suggest can be, and need to
be, obtained for a CMB-HD facility to achieve its lensing science goals [15].
We emphasize that a full demonstration that such foreground residuals can be
achieved using realistic simulations is a subject of ongoing research. We also
include baryon acoustic oscillation (BAO) data expected from the Dark Energy
Spectroscopic Instrument (DESI) [16], which serves to break a number of
parameter degeneracies.

Parameters are forecast with both a Fisher estimation code and a likelihood
plus Markov chain; we find that the resulting forecasts from both methods are
in agreement. We make publicly available the CMB-HD temperature, polarization,
and lensing ($TT,TE,EE,BB,\kappa\kappa$) signal and noise spectra, joint
covariance matrix, binning matrix, Fisher forecasting
software111https://github.com/CMB-HD/hdfisher, and
likelihood222https://github.com/CMB-HD/hdlike in order to enable further
forecasting. In addition, the CMB-HD likelihood has been integrated with the
latest version of Cobaya333https://cobaya.readthedocs.io [17] using
CAMB444https://camb.info/ [18, 19].

In Section II, we discuss the experimental configuration employed, our
generation of CMB and lensing noise spectra, and our method of delensing. In
Section III, we discuss the calculation of the joint covariance matrix. We
describe the Fisher estimation and likelihood plus Markov chain procedures for
forecasting cosmological parameters in Section IV. In Section V, we present
parameter forecasts for a CMB-HD survey and show the impact of foregrounds,
delensing, baryonic effects, and BAO. We discuss the implications of these
forecasts for discovering new light relic particles and probing inflation in
Section VI, and we conclude in Section VII.

##  II Generating Signal and Noise Spectra

Below we discuss our method of obtaining parameter forecasts for an ultra low-
noise, ultra-high-resolution CMB survey, such as CMB-HD [6, 7]. We focus in
particular on the $\Lambda$CDM model and minimal extensions. To assess the
additional constraining power of a CMB-HD experiment, we also provide
comparisons to experiments similar to the Simons Observatory (SO) [4] and
CMB-S4 [5]. We assume that the CMB datasets we forecast for will provide $TT$,
$TE$, $EE$, and $BB$ power spectra, as well as the $\kappa\kappa$ power
spectra from CMB lensing. In Section II.1, we detail the experimental
configurations assumed, and in Sections II.2 and II.3, we describe the
instrument noise and residual foreground models. In Section II.4, we describe
the computation of the CMB lensing power spectrum noise, and in Section II.5,
we explain the generation of the delensed CMB spectra.

###  II.1 Experimental Configuration

We create mock signal, noise, and covariance matrices for experiments similar
to CMB-HD, CMB-S4, SO, and DESI. Note that we do not aim to reproduce the
characteristics of these experiments exactly, but instead aim to give
approximate estimates based on the experimental configurations described below
and in Table 1.

For the CMB experiments we assume $60\%$ of the sky is observed. We also model
only the 90 and 150 GHz channels since they contain the bulk of the CMB
signal, and we assume the other frequencies will be utilized primarily for
foreground cleaning. We obtain the noise levels for HD-like, S4-like, and SO-
like surveys from [6, 5, 4], respectively, stressing again that we do not aim
to present official forecasts for each experiment, but rather to assess their
comparative constraining power. Throughout this work, “SO-like” and “S4-like”
surveys will refer to experimental configurations similar to the SO goal and
CMB-S4 wide-area surveys [4, 5]. We assume SO-like and S4-like experiments
will employ 5-meter dishes and a CMB-HD-like survey will use 30-meter dishes,
and we determine the resolution of each facility from these specifications.

Regarding the CMB multipole ranges, we assume in our forecasts that a CMB-HD-
like survey can measure multipoles in the range of $\ell\in[1000,\ 20000]$ for
$TT$, $TE$, $EE$. We assume CMB-HD will not measure $\ell<1000$ since these
multipoles are a challenge for a 30-meter dish555While the baseline CMB-HD
design does not assume multipoles below $\ell=1000$ will be measured by CMB-
HD, we note that there are interesting science cases that benefit from CMB-HD
measuring the $BB$ spectrum down to $\ell=100$, e.g. [20]; such low-$\ell$
measurements will be strived for if they can be achieved technologically.;
however, since both SO and CMB-S4 will measure the $TT$, $TE$, and $EE$
spectra to the sample-variance limit for $\ell\in[30,1000]$ over the same
$60\%$ of the sky as CMB-HD, we extend the CMB-HD multipole range down to
$\ell=30$ for $TT$, $TE$, and $EE$. For the $BB$ spectrum in the multipole
range of $30\leq\ell<1000$, we use the anticipated $BB$ noise from an Advanced
Simons Observatory (ASO) type survey; for ASO we assume the same resolution
and sky area as SO, but 3.5 and 3.8 $\mu$K-arcmin white noise in temperature
for 90 and 150 GHz, respectively.

Exp. | $f_{\mathrm{sky}}$ | Freq., | Noise, | FWHM, | Multipole Range
---|---|---|---|---|---
GHz |  $\mu$K-arcmin | arcmin | $\ell,L$
HD | 0.6 | 90 | 0.7 | 0.42 | TT, TE, EE: [30, 20000]
150 | 0.8 | 0.25 | BB: [1000, 20000]
|  |  |  $\kappa\kappa$: [30, 20000]
S4 | 0.6 | 90 | 2.0 | 2.2 | TT: [30, 3000]
150 | 2.0 | 1.4 | TE, EE, BB: [30, 5000]
|  |  |  $\kappa\kappa$: [30, 3000]
SO | 0.6 | 90 | 5.8 | 2.2 | TT: [30, 3000]
150 | 6.3 | 1.4 | TE, EE, BB: [30, 5000]
|  |  |  $\kappa\kappa$: [30, 3000]
Table 1: For each CMB experiment considered, we list the white noise level in
temperature and beam full-width at half-maximum (FWHM) for each frequency used
in the forecast, the sky fraction ($f_{\mathrm{sky}}$), and the multipole
ranges used for the power spectra. For the SO-like and S4-like forecasts, we
only forecast for the large-aperture telescopes, and use noise levels from [4,
5], respectively. For the CMB-HD-like forecasts, we use noise levels from [7],
and we assume CMB-HD only measures $\ell\geq 1000$. Since SO will measure
$TT,TE,$ and $EE$ to the sample variance limit for $\ell<1000$ over the CMB-HD
sky area, we extend the CMB-HD multipole range down to $\ell=30$ for these
spectra; for the $BB$ spectrum, we use an ASO-like $BB$ spectrum and noise in
the multipole range of $30\leq\ell<1000$ (see text for details). We emphasize
that we do not aim to present official forecasts for each experiment, but
rather to assess the comparative constraining power varying a number of
parameters.

Given the resolution of SO and CMB-S4, we limit the $TE$, $EE$, and $BB$
multipole range to $\ell\leq 5000$. While we do not explicitly include
foreground residuals in our SO-like and S4-like forecasts, we only consider
multipoles of $\ell\leq 3000$ for $TT$, assuming higher $TT$ multipoles will
be foreground dominated. Similarly, in light of foreground contamination and
bias, we only consider $\kappa\kappa$ for $L\leq 3000$ for SO-like and S4-like
experiments. In contrast, we include $L\in[30,20000]$ for an HD-like
experiment, and use the $\kappa\kappa$ covariance matrix from [15] for
$L\in[4800,20000]$, as we discuss in more detail below. We also include in the
HD-like CMB power spectra, the foreground residual levels required for CMB-HD
to achieve its lensing measurement target, as specified in [15].

To simulate the mock DESI dataset, we assume the data is taken over 14,000
square degrees (approximately 35% of the sky) for the baseline galaxy survey
and the bright galaxy survey, as specified in [16]. Our mock DESI data
consists of distance ratio measurements, $r_{\mathrm{s}}/d_{V}(z)$, obtained
from CAMB [18, 19], and a covariance matrix for those measurements, which we
calculate. Here $r_{\mathrm{s}}$ is the comoving sound horizon and $d_{V}$ is
a combined distance measurement, which we define in Section III.2. We
calculate the covariance matrix using the forecasted fractional errors on
$H(z)r_{\mathrm{s}}$ and $r_{\mathrm{s}}/d_{A}(z)$ provided in [16], as we
discuss in more detail in Section III.2.

###  II.2 CMB Instrumental Noise Spectra

A first step towards generating theoretical CMB and CMB lensing power spectra
is to compute the expected instrumental noise on the CMB power spectra. We add
to these noise spectra for $TT$ residual extragalactic foregrounds in the case
of CMB-HD, as described in Section II.3. These mock CMB power spectra are then
used to calculate the expected lensing reconstruction noise (Section II.4),
which is then used to calculate the theoretical delensed CMB power spectra
(Section II.5). The lensed or delensed CMB spectra and the lensing spectrum,
along with the CMB and lensing noise spectra, are then used in the analytic
calculation of the covariance matrix (Section III). This covariance matrix is
then used to forecast parameter constraints either via a Fisher matrix
(Section IV.1) or a likelihood and Markov Chain (Section IV.2).

The power spectra of the CMB temperature and polarization fields and the
(projected) lensing potential are determined by the cosmological model that
describes our Universe. We use a flat $\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum
m_{\nu}$ model described by eight cosmological parameters with fiducial values
from the Planck baseline cosmological parameter constraints [1], which we list
in Table 2. We use CAMB to generate the lensed, unlensed666These are used in
the analytic covariance matrix calculation described in Section III., and
delensed CMB power spectra $C_{\ell}^{XY}$ for $XY\in[TT,\ TE,\ EE,\ BB]$ and
the lensing power spectrum $C_{L}^{\kappa\kappa}$. We use the 2016 version of
HMcode [21, 22] to calculate the non-linear matter power spectrum for a cold
dark matter model, and increase the accuracy of CAMB beyond its default
settings as discussed in detail in Appendix A. We also use CAMB to calculate
the theoretical BAO signal, i.e. $r_{s}/d_{V}(z)$, using the same fiducial
cosmology and the default CAMB accuracy settings.

Parameter |  Fiducial | Step | Prior
---|---|---|---
size (%)
$\Omega_{\mathrm{b}}h^{2}$ |  $0.02237$ | 1 | $[0.005,\ 0.1]$
$\Omega_{\mathrm{c}}h^{2}$ |  $0.1200$ | 1 | $[0.001,\ 0.99]$
$\ln(10^{10}A_{\mathrm{s}})$ |  $3.044$ | 0.3777The step size of 0.3% $\ln(10^{10}A_{\mathrm{s}})$ corresponds to an approximate step size of 1% on $A_{\mathrm{s}}$. | $[2,\ 4]$
$n_{\mathrm{s}}$ |  $0.9649$ | 1 | $[0.8,\ 1.2]$
$\tau$ |  $0.0544$ | 5 | $0.054\pm 0.007$
$100\theta_{\mathrm{MC}}$ |  $1.04092$ | 1 | $[0.5,\ 10]$
$H_{0}$ [km s-1 Mpc-1] |  $67.36$ | 1 | $[20,\ 100]$
$N_{\mathrm{eff}}$ |  $3.046$ | 5 | $[0.05,\ 10]$
$\sum m_{\nu}$ [eV] |  $0.06$ | 10 | $[0,\ 5]$
Table 2: Listed in the first two columns are the cosmological parameters
varied in the Fisher or MCMC forecasts and their fiducial values based on the
Planck baseline cosmological parameter constraints [1]. The third column lists
the step size for each parameter, as a percentage of the fiducial value, used
when calculating the numerical derivatives of the power spectra with respect
to that parameter. The last column lists the priors used for the MCMC
analysis; the prior on the reionization optical depth, $\tau$, is also applied
to the Fisher forecasts. Note that we sample over $100\theta_{\mathrm{MC}}$
and obtain $H_{0}$ as a derived parameter in the MCMC analysis. All priors are
uniform, with the exception of $\tau$, where we use a Gaussian prior from
Planck [1].

We generate noise spectra for each CMB experiment listed in Table 1 at 90 and
150 GHz. Each frequency $f$ has a root-mean-square white noise level
$\Delta^{TT}_{f}$ in temperature and a beam full-width at half-maximum of
$\theta^{\mathrm{FWHM}}_{f}$ given in Table 1 . We assume that the noise in
polarization is given by
$\Delta^{EE}_{f}=\Delta^{BB}_{f}=\sqrt{2}\Delta^{TT}_{f}$, and set
$\Delta^{TE}_{f}=0$ assuming uncorrelated noise in temperature and
polarization maps. For a given frequency, we model the beam-deconvolved noise
power spectrum as

$N_{\ell_{f}}^{XY}=\left(\Delta^{XY}_{f}\right)^{2}\exp\left[\frac{\ell(\ell+1)\left(\theta_{f}^{\mathrm{FWHM}}\right)^{2}}{8\ln
2}\right],$ (1)

for $XY\in[TT,\ EE,\ BB]$.

To model a CMB-HD-type survey, we add the anticipated residual foreground
power spectra (see Section II.3) to the temperature noise spectrum at each
frequency, so that
${N_{\ell_{f}}^{TT}}^{\mathrm{HD}}=N_{\ell_{f}}^{TT}+C_{\ell_{f}}^{FG}$. Since
we assume that CMB-HD will not measure $\ell<1000$ (and does not need to for
$TT,TE,EE$ as discussed in Section II.1), we use an ASO-like $BB$ noise
spectrum for $\ell\in[30,1000)$. We assume polarized foregrounds can be
removed to levels below the instrument noise for all CMB experiments
considered here by making use of the additional frequency coverage outside of
90 and 150 GHz that each experiment will have.

We coadd the noise spectra for 90 and 150 GHz for each experiment using
inverse-noise weights, $W_{\ell_{f}}^{XY}=(N_{\ell_{f}}^{XY})^{-1}$, such that
the coadded noise power spectrum is given by

$N_{\ell}^{XY}=\frac{\sum_{f}\left(W_{\ell_{f}}^{XY}\right)^{2}N_{\ell_{f}}^{XY}}{\left(\sum_{f}W_{\ell_{f}}^{XY}\right)^{2}}.$
(2)

The total CMB power spectra are obtained by adding the coadded noise spectra
to the theoretical CMB power spectra, i.e.
$C_{\ell}^{XY,\mathrm{tot}}=C_{\ell}^{XY}+N_{\ell}^{XY}$. For parameter
forecasts (Sections IV.1 and IV.2), we use a binned signal spectra given by

$C_{\ell_{b}}^{XY}=\sum_{\ell}M_{b\ell}C_{\ell}^{XY},$ (3)

for $XY\in[TT,\ TE,\ EE,\ BB,\ \kappa\kappa]$, where $M_{b\ell}$ is the
binning matrix and $C_{\ell_{b}}^{XY}$ is the binned spectrum with bin centers
$\ell_{b}$. The binning is chosen to capture the CMB acoustic peaks on scales
of $\ell\lesssim 5000$ and uses a uniform bin width of $\Delta\ell=300$ on
smaller scales. The first four panels of Fig. 1 show the CMB theory,
instrumental noise, residual extragalactic foregrounds, and binned forecasted
error bars for CMB-HD, plotting the binned diagonals of the covariance matrix
(see Section III). In the last panel, we show the lensing theory power
spectrum, the lensing reconstruction noise (see Section II.4), and the binned
forecasted lensing error bars.

Figure 1:  Expected CMB $TT$, $TE$, $EE$, $BB$, and $\kappa\kappa$ power
spectra for a CMB-HD-type survey, with instrument noise and residual
extragalactic foregrounds coadded from $90$ and $150$ GHz. We show the theory
spectra (blue), the per-mode beam-deconvolved instrument noise (purple), the
residual temperature foregrounds (orange), and the expected CMB-HD error bars,
plotting the diagonals of the binned covariance matrix (red) (see Sections
II.1, II.2, and II.3 for details). We also show the noise curve for the
lensing power spectrum, $\kappa\kappa$, (pink), which is a combination of
noise from the primordial CMB, instrument, and residual foregrounds (see
Section II.4 for details). Note that in many cases the error bars are smaller
than the point indicating the bin center.

###  II.3 Residual Extragalactic Foregrounds

The observed CMB is contaminated by astrophysical foregrounds, which can be
both Galactic and extragalactic in origin. Galactic foregrounds dominate on
relatively large-scales, i.e. $\ell<500$, and are frequency dependent. As
mentioned above, we assume CMB-HD only measures multipoles above $\ell=1000$,
and will use lower multipoles from SO for $TT,TE$, and $EE$, which will
already be sample variance limited over that range. SO will have observations
over at least six different frequencies, and we assume in this work that SO
can remove the Galactic foreground contribution to below the CMB sample
variance limit for $TT,TE$, and $EE$. We assume CMB-HD will use the $BB$
measurement from ASO below $\ell=1000$, and we assume that ASO will remove
Galactic foregrounds to below their noise levels for $BB$ over this range.

For extragalactic foregrounds, the six times higher resolution and more than
three times lower noise of CMB-HD compared to precursor wide-area CMB surveys,
in addition to seven frequency channels, will allow for much greater
foreground subtraction. These extragalactic foregrounds include the thermal
and kinetic Sunyaev-Zel’dovich effects (tSZ and kSZ) from hot gas in galaxy
clusters, groups, and the intergalactic medium [23, 24, 25], the cosmic
infrared background (CIB) from dusty star-forming galaxies, and the radio
emission from active galactic nuclei. For CMB-HD, we assume the foreground
levels obtained in [15]; in particular, these residual foreground levels
assume a CMB-HD-type experiment will template-subtract CIB sources above 0.03
and 0.008 mJy at 150 and 90 GHz, respectively, in part exploiting CMB-HD
observations at 220 and 280 GHz to find sources and extrapolate their fluxes
to lower frequencies.888In the case of the CIB, see [15] for discussion of how
these flux levels were obtained, accounting for confusion from blended sources
and CIB spectral index uncertainties. These levels also assume radio sources
above 0.03 and 0.04 mJy at 150 and 90 GHz, respectively, are template
subtracted, as well as all tSZ sources detected at $3\sigma$ significance or
higher. We do not assume the kSZ effect from reionization can be removed,
since that is quite challenging given its frequency independence and high
redshift; however, we do assume the reionization kSZ only adds Gaussian noise
uncorrelated with the lensing signal. We assume that the late-time kSZ effect
can be removed from the CMB-HD maps through innovative “de-kSZ-ing” techniques
that are currently in development [26], and do not include it here. We stress
again that the residual foreground levels assumed are what preliminary studies
suggest is possible and necessary for CMB-HD to achieved an unbiased CMB
lensing power spectrum over the lensing multipole range of $L\in[5000,20000]$
[15]. However, we leave it to future work to robustly demonstrate with
realistic simulations and an end-to-end pipeline that these residual
foreground levels can be achieved in practice.

The total residual extragalactic foreground power from the sources described
above is shown as the dashed orange curve in the top panel of Fig. 1 for the
coadded 90 and 150 GHz temperature power spectra; this residual foreground
power is included in the total noise spectra for the CMB-HD temperature power
spectrum. We include only instrumental noise for the $TE,EE$, and $BB$ spectra
since we assume that the combination of high resolution and seven frequency
channels for CMB-HD will be sufficient to reduce polarized extragalactic
foregrounds to residual levels below the instrument noise.

As mentioned in Section II.1, while we do not explicitly include residual
extragalactic foreground levels in the temperature spectra for SO and CMB-S4,
we only include multipoles below $\ell=3000$ in those spectra, assuming higher
multipoles will be dominated by foreground uncertainty.

###  II.4 CMB Lensing Power Spectrum Noise

Lensing breaks the Gaussianity of the primordial CMB, introducing mode-
coupling in the CMB on different scales. A quadratic estimator [27, 28, 29,
30] takes advantage of this lensing-induced mode-coupling of the CMB to
reconstruct the lensing potential from pairs of CMB maps that are filtered to
isolate this mode-coupling. The lensing reconstruction can then be used to
estimate the lensing potential power spectrum, which is the connected (non-
Gaussian) part of the CMB four-point function. The disconnected (Gaussian)
part of the four-point function, called the Gaussian $N_{L}^{(0)}$ bias,
arises from the primordial CMB and instrument noise even in the absence of
lensing [31]; we assume that this will be subtracted with traditional
realization-dependent (RDN0) subtraction techniques [32].

We use two types of quadratic estimators in this work. The traditional
estimator from [27, 28, 29], henceforth called H&O, which uses the same CMB
multipole ranges in both CMB maps that go into the estimator. The other
quadratic estimator we use is from [30], henceforth called HDV, which allows
one map to contain only large scales ($\ell<2000$), called a gradient map, and
one map to contain only small scales ($\ell\in(5000,20000)$). The HDV
estimator is well suited to returning the small-scale signal, and allows one
to minimize bias from non-Gaussian astrophysical foregrounds by utilizing two
CMB maps that do not overlap in scales [15].

While we can remove the Gaussian noise bias by RDN0 subtraction techniques,
the Gaussian part of the four-point function still contributes noise to the
lensing power spectrum which can not be completely removed.999We note that
split-based lensing estimators can remove the contribution to $N_{L}$ from
instrument noise [33], but not from the primordial CMB. In this work, we do
not assume split-based estimators are used, but their use is a potential
avenue for gaining improved signal-to-noise. We hereafter call this lensing
spectrum noise $N_{L}^{\kappa\kappa}$ to distinguish it from the CMB noise
spectra $N_{\ell}^{XY}$ for $XY\in[TT,\ EE,\ BB]$.

We use the public CLASS delens package
[12]101010https://github.com/selimhotinli/class_delens to calculate
$N_{L}^{\kappa\kappa}$ for the H&O estimator [29].111111We cross checked that
the CLASS delens package returns the same $N_{L}^{\kappa\kappa}$ as the public
tempura package (https://github.com/simonsobs/tempura). We use the CLASS
delens package because it can iteratively delens all the estimator
combinations, not just $EB$, as we discuss in Section II.5. For SO-like and
CMB-S4-like experiments, we use the H&O estimator to generate
$N_{L}^{\kappa\kappa}$ curves using the CMB multipole ranges specified in
Table 1 and four-point lensing estimators from pairs of $TT,TE,EE,EB$, and
$TB$ maps, noting that the other possible four-point combinations contribute
less significantly to the overall signal-to-noise. For all the H&O
$N_{L}^{\kappa\kappa}$ curves, discussed here and below, we reduce their noise
levels further by iterative delensing, which we describe in Section II.5. We
coadd these $N_{L}^{\kappa\kappa}$ spectra to obtain a minimum variance (MV)
noise curve, $MV\leavevmode\nobreak\ N_{L}^{\kappa\kappa}$.

For the CMB-HD-like experiment, for lensing multipoles below $L=5000$, we
generate $N_{L}^{\kappa\kappa}$ curves using all CMB multipoles within
$\ell\in[30,20000]$ for $T,E$, and $B$121212Note that we use ASO-like noise
for the $BB$ spectrum for $\ell<1000$, and since CLASS delens does not allow
different noise levels for $E$ and $B$ for the $EB$ estimator, for $EB$ we set
both $E$ and $B$ noise levels to that of ASO below $\ell<1000$., but note that
there is minimal change to $N_{L}^{\kappa\kappa}$ for $L<5000$ if we limit the
$T$ spectra to $\ell\in[30,10000]$ due to foreground concerns. We calculate
$N_{L}^{\kappa\kappa}$ with the H&O estimator for $TT,TE,EE,EB$, and $TB$
maps, reducing these noise levels further by iterative delensing, and coadd
these to obtain an $MV\leavevmode\nobreak\ N_{L}^{\kappa\kappa}$. We assume
the $C_{L}^{\kappa\kappa}$ covariance matrix for $L\in[30,5000]$ has only
diagonal contributions since the lensing convergence is assumed to be Gaussian
on relatively large scales.

Figure 2: We show the theoretical lensing convergence power spectrum,
$\kappa\kappa$, (black solid) along with the CMB-HD minimum-variance (MV)
lensing power spectrum noise, $N_{L}^{\kappa\kappa}$, (red dashed). As
discussed in Section II.4, the latter is a combination of the estimators of
[28] (H&O) and [30] (HDV); in the lensing multipole range $L\in[30,5000]$, we
use the H&O MV reconstruction noise from the $TT$ (blue) and $TE$, $EE$, $EB$,
and $TB$ estimators (the MV of the latter four estimators are shown in
purple). For each of these estimators, we assume iterative delensing will be
done, and lower the noise curves accordingly following the method of [12] and
summarized in Section II.5. In the range $L\in[5000,20000]$, to reduce bias
from residual extragalactic foregrounds, we replace the $N_{L}^{\kappa\kappa}$
from the H&O $TT$ estimator with that from the HDV $TT$ estimator, which we
obtain from the diagonal elements of the simulation-based covariance matrix of
[15] (cyan solid); to this we coadd the H&O MV $N_{L}^{\kappa\kappa}$ from
$TE$, $EE$, $EB$, and $TB$. For comparison, we also show the analytically-
derived HDV $TT$ $N_{L}^{\kappa\kappa}$ (cyan dotted); the excess variance of
the simulation-based HDV $N_{L}^{\kappa\kappa}$ compared to the analytic one
is likely due to higher-order lensing corrections and non-Gaussian
fluctuations of the matter power spectrum that are significant on these small
scales and are captured by the simulations, as discussed in [34, 15]. See
Section II.4 for full details.

In Fig. 2, we show $N_{L}^{\kappa\kappa}$ for CMB-HD from the H&O estimator
from $TT$ alone (solid dark blue curve), and from everything other than $TT$
(i.e. minimum variance coadd of $TE,EE,EB,TB$) (solid purple curve). From the
comparison of these curves, we see that the polarization estimators dominate
the signal-to-noise below $L=2500$; however, on smaller scales, the $TT$
estimator has significantly lower noise than all the other estimators
combined. For this reason, on small scales, one must be concerned about the
impact of astrophysical foregrounds in the temperature maps and how they may
bias the lensing power spectrum. We also show in this figure the
$N_{L}^{\kappa\kappa}$ for the HDV quadratic estimator using only $TT$ and
$\ell\in[30,2000]$ in the gradient map and $\ell\in[5000,20000]$ in the small-
scale map (dotted cyan curve). This noise curve was made using the symlens
package131313https://github.com/simonsobs/symlens. We see that the HDV curve
matches the H&O curve well on small-scales; however, it is not as useful on
large scales. We note that all the $N_{L}^{\kappa\kappa}$ curves shown in Fig.
2 were calculated including both the instrument noise and residual foregrounds
described in Sections II.2 and II.3 for the CMB spectra.

For CMB-HD, on small lensing scales from $L\in[5000,20000]$, we use the
simulation-based lensing covariance matrix from [15] for $TT$, which
incorporates the HDV quadratic estimator with the $\ell$ ranges for large-
scale and small-scale CMB maps given above, as well as higher-order lensing
corrections and non-Gaussian fluctuations of the matter power spectrum that
are significant on these small scales and are captured by the simulations. The
solid cyan curve in Fig. 2 shows the lensing noise obtained from the diagonals
of this covariance matrix. This $C_{L}^{\kappa\kappa}$ covariance matrix also
includes off-diagonal contributions (which we show in Fig. 6); in generating
this covariance matrix, simulation-based calculations of the RDN0 and $N_{1}$
biases (the latter is a higher-order correction) are subtracted, which reduces
the off-diagonal correlations significantly, as discussed in [15] and [34].
Restricting the HDV gradient map to $\ell<2000$ is also important for removing
the $N_{2}$ bias [34, 35]. For the final CMB-HD $N_{L}^{\kappa\kappa}$ spectra
in the range of $L\in[5000,20000]$, we coadd the simulation-based $TT$ HDV
$N_{L}^{\kappa\kappa}$ with the H&O $N_{L}^{\kappa\kappa}$ from $TE,EE,EB,$
and $TB$. For $L\in[30,5000]$, we use the H&O minimum variance
$N_{L}^{\kappa\kappa}$ coadding $TT,TE,EE,EB,$ and $TB$ noise curves. We show
this CMB-HD MV $N_{L}^{\kappa\kappa}$ curve in Fig. 1 (dashed pink) and in
Fig. 2 (dashed red).

For CMB-HD, the $EB$ estimator dominates below $L=2500$, and foreground biases
are a negligible concern. Above $L=5000$, we attempt to be conservative by
adopting the HDV estimator and a simulation-based covariance matrix; we note
that there are other more optimal estimators that can be used to reconstruct
lensing at small scales that may yield a higher signal-to-noise ratio and
which also hold promise for mitigating foreground bias [36, 37, 9, 38]. In the
regime of $L\in[2500,5000]$, where the $TT$ estimator dominates, we are
optimistic that foreground bias can be mitigated, but emphasize that this
still needs to be demonstrated with realistic simulations; as mentioned above,
the $TT$ spectrum can be cut to only include $\ell\in[30,10000]$ as opposed to
$\ell\in[30,20000]$, with no loss of lensing signal-to-noise for $L<5000$.

###  II.5 Creation of Delensed CMB Spectra

Figure 3: The residual lensing power, $C_{L}^{\kappa\kappa,\mathrm{res}}$,
remaining after delensing (Eq. 4), shown as a fraction of the full lensing
power spectrum, $C_{L}^{\kappa\kappa}$, for SO, CMB-S4, and CMB-HD-like
experiments with the specifications listed in Table 1. For CMB-HD, the change
near $L=5000$ is due to using different lensing estimators above and below
that scale, as discussed in Section II.4. Below $L=650$, CMB-HD removes over
90% of the lensing power.

Given the ultra low noise of a CMB-HD-type survey, it will be possible to
remove a significant amount of the lensing in the CMB spectra, a procedure
called delensing. Delensing has the effect of tightening cosmological
parameter constraints because it sharpens the acoustic peaks, removes
uncertainty in the lensing realization that distorted the primordial CMB, and
reduces lensing-induced off-diagonal power spectrum covariances [11, 12]. In
order to achieve tighter parameter constraints via delensing, it is important
to replace the lost lensing information in the CMB power spectra with that
from the reconstructed lensing power spectrum by combining them in the
likelihood, properly accounting for the delensed CMB and CMB lensing power
spectrum covariances [11].

Figure 4: Shown is the amount of lensing power remaining in CMB-HD $TT$, $TE$,
and $EE$ power spectra after delensing, $\Delta C_{\ell}^{\mathrm{delensed}}$,
as a fraction of the total amount of lensing present in the CMB spectra,
$\Delta C_{\ell}^{\mathrm{lensed}}$. $\Delta C_{\ell}$ is defined as the
difference between the lensed/delensed CMB power spectra and the unlensed
power spectra. Since both $\Delta C_{\ell}$ curves cross zero at the same
multipoles, we bin them before taking their ratios, and show the bin centers
as points. The lensing in the CMB-HD spectra is reduced on all scales, and
about 90% of the lensing is removed in the CMB spectra for $\ell<2000$
(indicated by points below the dotted line). See Section II.5 for details.

To predict the amount of lensing that can be removed from the CMB power
spectra for each experiment discussed in Section II.1, we first calculate the
residual lensing power that will remain in the CMB maps after delensing. We
then lens unlensed CMB spectra by this remaining lensing power using CAMB,
following [8]. The residual lensing power,
$C_{L}^{\kappa\kappa,\mathrm{res}}$, is calculated by subtracting the Wiener-
filtered lensing power spectrum from the total theoretical
$C_{L}^{\kappa\kappa}$, as follows

$C_{L}^{\kappa\kappa,\mathrm{res}}=C_{L}^{\kappa\kappa}\left(1-\frac{C_{L}^{\kappa\kappa}}{C_{L}^{\kappa\kappa}+N_{L}^{\kappa\kappa}}\right),$
(4)

where $N_{L}^{\kappa\kappa}$ is the expected lensing reconstruction noise for
the MV estimator described in Section II.4. This removes lensing on scales
where the lensing power spectrum is detected with high signal-to-noise ratio,
but not on scales where the lensing reconstruction is too noisy to estimate
and remove the lensing signal. We show the residual lensing power as a
fraction of the total power for SO, CMB-S4, and CMB-HD-like experiments in
Fig. 3. We see that below $L=650$, delensing with CMB-HD removes over 90% of
the lensing power. The abrupt change near $L=5000$ is due to using different
lensing estimators above and below that scale, as discussed in Section II.4.

As mentioned in Section II.4, for all of the H&O $N_{L}^{\kappa\kappa}$ curves
discussed above, we reduce their noise levels further by iterative delensing.
Iterative delensing of the $EB$ lensing estimator is traditionally assumed in
forecasts of future CMB experiments [5]. The reason is that, assuming no
primordial $B$-mode power (i.e., assuming the $B$-mode power is entirely due
to lensing), there are no unlensed $B$-modes to contribute to the variance of
the lensing reconstruction. Thus the $EB$ lensing reconstruction noise can be
reduced by delensing, which lowers the $B$-mode power. Then the delensed $E$
and $B$ maps can be used to obtain a reconstruction of the residual lensing,
which can be used to further delens the $E$ and $B$ maps; this procedure can
be iterated on until converged [39, 40, 41]. On small scales, where most of
the CMB power in the temperature and $E$-mode polarization maps is due to
lensing, the amount of delensing can be improved by a similar iterative
approach [12], where all the estimators are iterated instead of only $EB$.

Figure 5: We show the amount of lensing that can be removed from any $BB$
power spectrum via delensing with CMB-HD. In particular, we show the lensed
and delensed CMB $BB$ power spectra (top panel), and their ratio (lower
panel), in the region where delensing is most effective ($\ell<2000$). The
bottom panel shows that CMB-HD delensing can remove over 90% of the lensing
signal in the $BB$ spectrum for $\ell<1000$; this is the case even when
assuming CMB-HD will not measure $\ell<1000$ and instead using ASO-like $BB$
spectra and instrument noise for $\ell\in[30,1000)$, as discussed in Section
II.2. This is because the higher CMB multipoles of CMB-HD are being used to
reconstruct low $L$ lensing multipoles used for delensing.

As mentioned in Section II.4, we use the CLASS delens code [12] to calculate
the expected lensing reconstruction noise when using the iterative delensing
procedure, lowering the final $MV$ lensing noise spectrum, as well as
increasing the amount of delensing of the CMB spectra, accordingly. We show
the effect of iterative delensing on the CMB $TT$, $TE$, and $EE$ power
spectra in Fig. 4, and on the CMB $BB$ power spectrum in Fig. 5. We see that
CMB-HD removes about 90% of the lensing in the $TT$, $TE$, and $EE$ spectra
for $\ell<2000$, and 90% of the lensing in the $BB$ spectra for $\ell<1000$.

In this work we assume that the bias that arises from internally delensing the
acoustic peaks [42, 10] can either be removed following a method similar to
[8] or will not arise as is the case for forward modeling methods [9]. We
leave for future work the demonstration of the mitigation of this bias.

##  III Covariance Matrix Calculation

In addition to the signal and noise spectra described in Section II, a full
covariance matrix of the expected data is required to forecast cosmological
parameter constraints.

###  III.1 CMB Covariance Matrix

To construct the covariance matrix for the CMB experiments discussed above, we
assume we will have $TT,TE,EE,BB$, and $\kappa\kappa$ spectra for each
experiment for the multipole ranges listed in Table 1. We analytically
calculate the covariance matrix following [12] as described below.

For the diagonal components of the covariance matrix, for the blocks giving
the variance between different CMB spectra (see Fig. 6), we calculate

$\mathbb{C}_{\ell_{1}\ell_{2},\mathrm{G}}^{XY,WZ}=\frac{\delta_{\ell_{1}\ell_{2}}f_{\mathrm{sky}}^{-1}}{2\ell_{1}+1}\left[C_{\ell_{1}}^{XW,\mathrm{tot}}C_{\ell_{1}}^{YZ,\mathrm{tot}}+C_{\ell_{1}}^{XZ,\mathrm{tot}}C_{\ell_{1}}^{YW,\mathrm{tot}}\right],$
(5)

Here $C_{\ell}^{\mathrm{tot}}$ is the total CMB signal plus noise spectrum,
(i.e. $C_{\ell}^{\mathrm{tot}}=C_{\ell}+N_{\ell})$, and combinations of any
two of $X$, $Y$, $Z$, and $W$ can each be $TT,TE,EE$, or $BB$. The sky
fraction observed by each experiment is indicated by $f_{\mathrm{sky}}$. We
label this term given by Eq. 5 with the subscript G for “Gaussian” to indicate
that it arises from the Gaussian fluctuations of the lensed or delensed CMB
field.

As discussed in Section II.4, lensing induces mode coupling between different
multipoles of the unlensed CMB field, giving the originally Gaussian field a
non-Gaussian structure. The uncertainty in the amount of mode coupling,
arising from uncertainty in either the lensing potential or the unlensed CMB
field due to sample variance, generates off-diagonal elements in the CMB x CMB
blocks [43]. Following [12], we calculate these off-diagonal terms as

$\begin{split}\mathbb{C}_{\ell_{1}\ell_{2},\mathrm{NG}}^{XY,WZ}&=\sum_{L}\frac{\partial
C_{\ell_{1}}^{XY}}{\partial
C_{L}^{\phi\phi}}\frac{2f_{\mathrm{sky}}^{-1}}{2L+1}\left(C_{L}^{\phi\phi}\right)^{2}\frac{\partial
C_{\ell_{2}}^{WZ}}{\partial C_{L}^{\phi\phi}}\\\
&\quad+\sum_{\ell}\left[\frac{\partial C_{\ell_{1}}^{XY}}{\partial
C_{\ell}^{XY,\mathrm{u}}}\mathbb{C}_{\ell\ell,\mathrm{G}}^{XY,WZ;\mathrm{u}}\frac{\partial
C_{\ell_{2}}^{WZ}}{\partial
C_{\ell}^{WZ,\mathrm{u}}}\right](1-\delta_{\ell_{1}\ell_{2}}),\end{split}$ (6)

where

$\mathbb{C}_{\ell\ell,\mathrm{G}}^{XY,WZ;\mathrm{u}}=f_{\mathrm{sky}}^{-1}\frac{\left(C_{\ell}^{XW,\mathrm{u}}C_{\ell}^{YZ,\mathrm{u}}+C_{\ell}^{XZ,\mathrm{u}}C_{\ell}^{YW,\mathrm{u}}\right)}{2\ell+1}$
(7)

is the sample variance of the unlensed CMB spectra $C_{\ell}^{\mathrm{u}}$,
and again combinations of any two of $X$, $Y$, $Z$, and $W$ can each be
$TT,TE,EE$, or $BB$.141414Note that when a $BB$ combination is involved, we
substitute the unlensed $EE$ spectrum in the second term of Eq. 6 and in Eq.
7, following [12, 44]. Similarly
$2f_{\mathrm{sky}}^{-1}(C_{L}^{\phi\phi})^{2}/(2L+1)$ is the sample variance
of the lensing potential spectrum. The lensing potential power spectrum,
$C_{L}^{\phi\phi}$, is related to the lensing deflection power spectrum,
$C_{L}^{dd}$, and the lensing convergence power spectrum,
$C_{L}^{\kappa\kappa}$, by

$C_{L}^{\kappa\kappa}=\frac{L(L+1)}{4}C_{L}^{dd}=\frac{L^{2}(L+1)^{2}}{4}C_{L}^{\phi\phi}.$
(8)

The subscript NG in Eq. 6 indicates that it is the “non-Gaussian” component.
The full covariance matrix for the CMB x CMB blocks is given by

$\mathbb{C}_{\ell_{1}\ell_{2}}^{XY,WZ}=\mathbb{C}_{\ell_{1}\ell_{2},\mathrm{G}}^{XY,WZ}+\mathbb{C}_{\ell_{1}\ell_{2},\mathrm{NG}}^{XY,WZ}.$
(9)

For the CMB lensing x CMB lensing block of the covariance matrix, the variance
of the lensing potential spectrum, after RDN0 subtraction, is approximately
Gaussian [44] and is given by

$\mathbb{C}_{L_{1}L_{2}}^{\phi\phi,\phi\phi}=\delta_{L_{1}L_{2}}\frac{2f_{\mathrm{sky}}^{-1}}{2L_{1}+1}\left(C_{L_{1}}^{\phi\phi}+N_{L_{1}}^{\phi\phi}\right)^{2}.$
(10)

For CMB-HD, we use Eq. 10 for $L<5000$, and for $L\geq 5000$ we replace the
analytic covariance matrix with a simulation-based one from [15] when
including foregrounds, or from [34] when not including foregrounds, both of
which use the HDV $TT$ estimator as discussed in Section II.4. Note
importantly that the simulation-based covariance matrices have off-diagonal
terms, as discussed in more detail in Section II.4, and as shown in Fig. 6 for
the $\kappa\kappa$ block.151515Note that on large-scales the lensing field is
approximately Gaussian, with only diagonal contributions to the covariance
matrix. In addition, RDN0 subtraction removes most residual off-diagonal
contributions [44]. However, a more accurate treatment would use a simulation-
based covariance matrix for all scales.

Figure 6: Left: We show the binned analytic correlation matrices for a CMB-HD
survey between the lensed CMB $TT$, $TE$, $EE$, and $BB$ power spectra (first
four rows), and between the lensed CMB and lensing $\kappa\kappa$ power
spectrum (last row). We set the diagonals of each block to zero to highlight
the off-diagonal elements, which arise due to the non-Gaussian covariances, as
discussed in Section III. We also show the simulation-based
$\kappa\kappa\times\kappa\kappa$ covariance matrix from [15] which we use on
scales $L,L^{\prime}\gtrsim 5000$; note that this is shown with a different
color scale than the analytic correlation matrices. Right: We show the
corresponding binned correlation matrices for the delensed CMB-HD spectra. We
find that the correlation coefficients for the off-diagonal terms of the
binned analytic correlation matrices (not including
$\kappa\kappa\times\kappa\kappa$) are less than $\pm 1\%$ everywhere.

We calculate the correlations between the lensing potential spectra and the
CMB spectra from Eq. 36 in [44], which can also be obtained by extending Eq. 5
in [12] to allow either $XY$ or $WZ$ to be $\phi\phi$, as

$\mathbb{C}_{L_{1}\ell_{2}}^{\phi\phi,XY}=\frac{2}{2L_{1}+1}\left(C_{L_{1}}^{\phi\phi}\right)^{2}\frac{\partial
C_{\ell_{2}}^{XY}}{\partial C_{L_{1}}^{\phi\phi}}.$ (11)

This covariance arises because the same lenses that make up the lensing
potential are also responsible for lensing the CMB spectra.

The terms given by Eqs. 5, 6, 10, and 11 are the most important terms in the
analytic covariance matrix calculation, with the diagonal/Gaussian terms
having the largest impact on parameters. However, there are some additional
terms that are not included here [45, 46, 44], which we discuss below. The
first term we neglect is the super-sample covariance [46] arising from lenses
on scales larger than the survey area. Since all three experiments listed in
Table 1 will survey half the sky, this effect is expected to be negligible as
suggested by [46].

For the covariance between the lensing potential spectrum and the CMB spectra,
there are four terms given in Eq. 33 of [44]. As discussed in [44], the first
two terms cancel if one does RDN0 subtraction, which we assume will be done
for all the CMB experiments considered here, as mentioned in Section II.4. The
third term in Eq. 33 of [44] is Eq. 11 above, and we found that the fourth
term (Type B Primary Trispectrum) had negligible impact on the parameter
forecasts and is thus omitted.

In this work, we neglect a first-order correction to the lensing power
spectrum known as the “N1 bias” [47, 35]. While this term is called a “bias”,
it is in fact a higher-order correction to the lensing signal that can be
calculated from first principles similar to the zeroth-order lensing spectrum.
This N1 signal is of similar or larger amplitude than the standard CMB power
spectrum at high multipoles, and its inclusion may yield a higher lensing
signal-to-noise ratio than we forecast here. We leave further exploration of
this to future work, and note that ignoring this term is likely conservative.

We use the CLASS delens package [12] to evaluate the derivatives $\partial
C_{\ell}^{XY}/\partial C_{L}^{\phi\phi}$ and $\partial C_{\ell}^{XY}/\partial
C_{\ell^{\prime}}^{XY,\mathrm{u}}$ in Eq. 6 and 11.161616We increase the CLASS
precision parameters beyond the default values to the following:
accurate_lensing = 1, delta_l_max = 2000, perturbations_sampling_stepsize =
0.05, k_max_tau0_over_l_max = 50, and P_k_max_1/Mpc = 500 to be consistent
with the accuracy settings used for CAMB (see Appendix A.) With these
derivatives, we assemble the final covariance matrix combining the terms in
Eqs. 5, 6, 10, and 11. We cross checked the CLASS delens $\partial
C_{\ell}^{XY}/\partial C_{L}^{\phi\phi}$ derivatives by comparing to those we
obtained using the lenscov171717https://github.com/JulienPeloton/lenscov
package [44] up to $\ell,L=10,000$. Since lenscov does not calculate
covariances involving delensed CMB spectra, nor does it compute the
derivatives $\partial C_{\ell}^{XY}/\partial
C_{\ell^{\prime}}^{XY,\mathrm{u}}$, we constructed covariance matrices from
the CLASS delens and lenscov packages using lensed CMB spectra and Eqs. 5, 6,
10, and 11 (excluding the second term in Eq. 6). We confirmed that the
forecasted parameter errors obtained from either covariance matrix were
consistent (within $<0.05$% for each parameter).

We generated covariance matrices for both lensed and delensed CMB spectra, by
setting the $C_{\ell_{1}}$ in Eq. 5 to either lensed or delensed spectra, with
the latter obtained as described in Section II.5. For the derivatives in Eqs.
6 and 11, CLASS delens provides an option to output them for delensed CMB
spectra, given the lensing spectrum noise, $N_{L}^{\kappa\kappa}$.

We bin our final covariance matrix using the same binning matrix, $M_{b\ell}$,
as used for the one-dimensional spectra (see Eq. 3 of Section II.2), so that
we obtain a binned covariance matrix given by

$\mathbb{C}_{\ell_{b}\ell_{b^{\prime}}}^{XY,WZ}=\sum_{\ell\ell^{\prime}}M_{b\ell}\mathbb{C}_{\ell\ell^{\prime}}^{XY,WZ}M_{b^{\prime}\ell^{\prime}}^{T},$
(12)

for $XY$ and $WZ$ in $[TT,\ TE,\ EE,\ BB,\ \kappa\kappa]$. We then form a
single block covariance matrix from the 25 binned blocks. We replace the
binned CMB-HD $\kappa\kappa\times\kappa\kappa$ covariance matrix on scales
$L_{1},\ L_{2}\gtrsim 5000$, with the simulation-based one from [15] when
including foregrounds, or [34] when not including foregrounds, as mentioned
above.

We show in Fig. 6 the binned analytic correlation matrices corresponding to
the lensed and delensed covariance matrices, calculated using Eq. 5, 6, 10,
and 11, for the CMB-HD experiment. We also show the binned simulation-based
$\kappa\kappa\times\kappa\kappa$ correlation matrix from [15] which is used on
scales $L_{1},\ L_{2}\gtrsim 5000$. We calculate the unbinned elements of the
correlation matrix, $\rho_{\ell_{1}\ell_{2}}^{XY,WZ}$, between
$C_{\ell_{1}}^{XY}$ and $C_{\ell_{2}}^{WZ}$ by [44, 12]

$\rho_{\ell_{1}\ell_{2}}^{XY,WZ}=\frac{\mathbb{C}_{\ell_{1}\ell_{2}}^{XY,WZ}}{\sqrt{\mathbb{C}_{\ell_{1}\ell_{1}}^{XY,XY}\mathbb{C}_{\ell_{2}\ell_{2}}^{WZ,WZ}}},$
(13)

where $XY$ and $WZ$ are $TT$, $TE$, $EE$, $BB$, or $\kappa\kappa$ (we have
used lowercase $\ell$’s for simplicity). We then bin this in the same way as
Eq. 12 to generate Fig. 6. For this figure, in order to highlight the off-
diagonal elements, the diagonal elements of each block were set to zero. In
the CMB x CMB blocks for lensed spectra at low multipoles (Fig. 6, left), we
see the familiar “checkerboard” pattern from the lensing-induced peak-
smoothing, as noted by [43, 45, 44, 11]. For the delensed correlation matrix
(Fig. 6, right), correlations in these blocks are reduced at low multipoles,
as also seen in [11, 12]; in fact they are nearly zero on scales where
delensing removes much of the lensing from the CMB spectra (as shown in Fig.
4). In the CMB x $\kappa\kappa$ blocks, correlations are visible between
large-scale (low-$L$) lensing and lensed CMB spectra for a wide range of
scales (up to $\ell\sim 5000$), as also discussed in [11]. The $\kappa\kappa$
x $\kappa\kappa$ blocks show off-diagonal terms from higher-order lensing
corrections and non-Gaussian fluctuations of the matter power spectrum that
are captured by simulations [15]. We find that in both the lensed and delensed
binned analytic correlation matrices (not including
$\kappa\kappa\times\kappa\kappa$), the correlation coefficients,
$\rho_{\ell_{1}\ell_{2}}^{XY,WZ}$, for the off-diagonal terms are less than
$\pm 1\%$.

###  III.2 BAO Covariance Matrix

We also construct a covariance matrix for BAO data from a DESI-like survey.
The BAO data measures the transverse and radial scales, $r_{s}/d_{A}(z)$ and
$H(z)r_{s}$, respectively, where $r_{s}$ is the comoving sound horizon at the
end of the baryon drag epoch (i.e. $r_{s}(z_{d})$), $d_{A}(z)$ is the angular
diameter distance to redshift $z$, and $H(z)$ is the expansion rate at
redshift $z$ [16]. These are combined into a distance measurement defined by
[48]

$d_{V}(z)\equiv\left[(1+z)^{2}d_{A}^{2}(z)\frac{cz}{H(z)}\right]^{1/3}.$ (14)

We assume a Gaussian (diagonal) covariance matrix for the DESI BAO data. We
calculate the uncertainty $\sigma_{j}$ on $r_{s}/d_{V}(z_{j})$ following the
approach of [49] as

$\begin{split}\sigma_{j}&\equiv\sigma\left(\frac{r_{s}}{d_{V}(z_{j})}\right)\\\
&=\frac{1}{3}\frac{r_{s}}{d_{V}(z_{j})}\sqrt{\left[\frac{\sigma\left(H(z_{j})r_{s}\right)}{H(z_{j})r_{s}}\right]^{2}+\left[2\frac{\sigma\left(d_{A}(z_{j})/r_{s}\right)}{d_{A}(z_{j})/r_{s}}\right]^{2}}.\end{split}$
(15)

We use CAMB to calculate the theoretical BAO signal $r_{s}/d_{V}(z_{j})$, at
each redshift $z_{j}$, for the Planck fiducial cosmology used in this work. We
use the forecasted fractional uncertainties on the quantities $H(z_{j})r_{s}$
and $d_{A}(z_{j})/r_{s}$, at redshifts $z_{j}$, given by Table 2.3 (for
$z\in[0.65,1.85]$) and Table 2.5 (for $z\in[0.05,0.45]$) in [16] for the DESI
experiment covering 14,000 square degrees of sky (i.e., the terms under the
square root in the second line of Eq. 15). The diagonal elements of the
covariance matrix are $\sigma_{j}^{2}$, with all off-diagonal elements set to
zero.

##  IV Methods for Parameter Forecasts

We use both a Fisher matrix method (described in Section IV.1) and a
likelihood plus Markov chain Monte Carlo (MCMC) method (described in Section
IV.2) to forecast cosmological parameter constraints from the combination of
CMB and BAO data. These forecasts are presented in Section V.

We forecast uncertainties on the six parameters of the $\Lambda$CDM model: the
physical baryon and cold dark matter densities, $\Omega_{\mathrm{b}}h^{2}$ and
$\Omega_{\mathrm{c}}h^{2}$, respectively; the primordial comoving curvature
power spectrum amplitude, $\ln(10^{10}A_{\mathrm{s}})$, and scalar spectral
index, $n_{\mathrm{s}}$, defined at the pivot scale $k_{0}=0.05$ Mpc-1; the
reionization optical depth, $\tau$; and the Hubble constant $H_{0}$. For MCMC
runs, we treat $H_{0}$ as a derived parameter, and instead sample over the
parameter $100\theta_{\mathrm{MC}}$, where $\theta_{\mathrm{MC}}$ is the
cosmoMC181818https://cosmologist.info/cosmomc/readme.html approximation to the
angular scale of the sound horizon at last scattering. We additionally vary
the effective number of relativistic species $N_{\mathrm{eff}}$ and the sum of
the neutrino masses $\sum m_{\nu}$. We refer to these eight parameters (listed
in Table 2) as the $\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$ model. Since
the CMB experiments considered here do not observe very large scales
($\ell<30$), we apply a prior on the optical depth $\tau$ of
$\sigma(\tau)=0.007$ from Planck [1].

Both the Fisher and MCMC methods assume that the CMB and BAO data can be
described by a Gaussian likelihood function
$\mathcal{L}(\hat{d}|\vec{\theta})$ for the probability distribution of the
data $\hat{d}$ given a set of model parameters $\vec{\theta}$. For the CMB
data described in Sections II and III.1, our likelihood function is given by

$-2\ln\mathcal{L}_{\mathrm{CMB}}\left(\hat{C}_{\ell_{b}}|\vec{\theta}\right)=\sum_{\ell_{b}\ell_{b^{\prime}}}\Delta
C_{\ell_{b}}(\vec{\theta})\mathbb{C}^{-1}_{\ell_{b}\ell_{b^{\prime}}}\Delta
C_{\ell_{b^{\prime}}}(\vec{\theta}),$ (16)

where $\Delta
C_{\ell_{b}}(\vec{\theta})=\hat{C}_{\ell_{b}}-C_{\ell_{b}}(\vec{\theta})$,
$\hat{C}_{\ell_{b}}$ is the binned (lensed or delensed) data spectra191919We
use theory spectra at our fiducial cosmology, with no scatter, to simulate the
data spectra. We verified that this yields the same error bars in our MCMC
analysis as we would get if we had added scatter to the spectra drawn from our
covariance matrix. with bin centers $\ell_{b}$, $C_{\ell_{b}}(\vec{\theta})$
is the binned theory spectra evaluated at a given set of cosmological
parameters $\vec{\theta}$, and $\mathbb{C}_{\ell_{b}\ell_{b^{\prime}}}$ is the
binned covariance matrix of the data. $\hat{C}_{\ell_{b}}$ and
$C_{\ell_{b}}(\vec{\theta})$ hold the binned $TT$, $TE$, $EE$, $BB$, and
$\kappa\kappa$ spectra (in the same format and order as the covariance matrix,
$\mathbb{C}_{\ell_{b}\ell_{b^{\prime}}}$).202020While the $BB$ power spectrum
does not add significantly to the parameter constraints, we include it for
completeness.

For the BAO data, given by $f_{j}\equiv r_{s}/d_{V}(z_{j})$ and described in
Section III.2, we use the likelihood function

$-2\ln\mathcal{L}_{\mathrm{BAO}}\left(\hat{f}_{j}\
|\vec{\theta}\right)=\sum_{j}\frac{\left[\hat{f}_{j}-f_{j}(\vec{\theta})\right]^{2}}{\sigma_{j}^{2}},$
(17)

where $\hat{f}_{j}$ is the BAO data with variance $\sigma_{j}^{2}$ at redshift
$z_{j}$, and $f_{j}(\vec{\theta})$ is the theory evaluated at the set of
cosmological parameters $\vec{\theta}$ [50].

We make both our Fisher and likelihood code publicly available, as mentioned
above, and describe them in more detail below.

###  IV.1 Fisher Matrix Estimation

Fisher matrices can be used as an alternative to MCMC runs for parameter
forecasts. While the Fisher approach is much quicker than a full MCMC run, it
assumes that the parameters follow Gaussian posterior distributions; this is a
good approximation for the six $\Lambda$CDM parameters, but may not be
suitable when allowing other parameters (such as neutrino mass) to vary. Thus
we cross check our Fisher matrix results with a likelihood/MCMC analysis and
find good agreement (see Appendix B).

To form a Fisher matrix, we calculate the change in the theoretical signal as
each parameter is varied (i.e. the numerical derivatives of spectra with
respect to each parameter); the covariance matrix of the data tells us how
well we can detect this change, and therefore how well we can constrain
parameters. We assume that the data can be described by a Gaussian likelihood,
such as the ones in Eqs. 16 and 17. We also assume that the likelihood is
maximized when evaluated at a fiducial set of parameters $\vec{\theta}_{0}$,
and Taylor expand about this point. The elements $F_{\alpha\beta}$ of the
Fisher matrix $F$ are given by the coefficients of the quadratic terms, i.e.

$F_{\alpha\beta}\equiv\left.-\left\langle\frac{\partial^{2}\ln\mathcal{L}}{\partial\theta_{\alpha}\partial\theta_{\beta}}\right\rangle\right|_{\vec{\theta}_{0}},$

where the indices $\alpha$ and $\beta$ correspond to specific parameters. This
characterizes how rapidly the likelihood changes as the parameters move away
from their fiducial values. The covariance matrix for the parameters is found
by taking the inverse of the Fisher matrix, i.e.
$\mathrm{cov}\left(\theta_{\alpha},\leavevmode\nobreak\
\theta_{\beta}\right)=\left(F^{-1}\right)_{\alpha\beta}$. The forecasted
marginalized error $\sigma_{\alpha}$ on parameter $\theta_{\alpha}$ is then
found from, $\sigma_{\alpha}=\sqrt{\left(F^{-1}\right)_{\alpha\alpha}}$.

For the CMB likelihood function given by Eq. 16, the elements of the Fisher
matrix are

$F_{\alpha\beta}^{\mathrm{CMB}}=\sum_{\ell_{b},\ell_{b}^{\prime}}\left.\left(\frac{\partial
C_{\ell_{b}}}{\partial\theta_{\alpha}}\mathbb{C}^{-1}_{\ell_{b}\ell_{b}^{\prime}}\frac{\partial
C_{\ell_{b}^{\prime}}}{\partial\theta_{\beta}}\right)\right|_{\vec{\theta}_{0}}.$
(18)

We use Eq. 18 to compute the CMB Fisher matrix for a given set of cosmological
parameters, numerically evaluating the derivatives with respect to parameters
using a finite difference method with the parameter step sizes listed in Table
2. We apply the Gaussian prior on $\tau$ by adding its inverse variance
$1/\sigma^{2}(\tau)$ to the corresponding element $F_{\tau\tau}$ of the Fisher
matrix [51].

For the BAO likelihood function given by Eq. 17, the elements of the Fisher
matrix are given by

$F_{\alpha\beta}^{\mathrm{BAO}}=\sum_{j}\left.\left(\frac{\partial
f_{j}}{\partial\theta_{\alpha}}\frac{1}{\sigma^{2}_{j}}\frac{\partial
f_{j}}{\partial\theta_{\beta}}\right)\right|_{\vec{\theta}_{0}},$ (19)

where the sum is taken over the redshifts $z_{j}$, and the derivatives
$\partial f_{j}/\partial\theta_{\alpha}$ are calculated in the same way as in
Eq. 18. To forecast parameter errors from the combination of CMB and BAO data,
we take the sum of the two Fisher matrices (since they are from independent
data sets), and then invert this sum to obtain the parameter covariance matrix
and marginalized error bars [51].

The Fisher matrix method described above gives the statistical uncertainties
on a set of parameters, assuming that the fiducial model used in the
calculation of the numerical derivatives accurately describes the true
universe. We can extend this method to calculate the expected bias to the
estimated parameter values due to an incorrect fiducial model, following the
approach described in [52, 53, 54] and adopted by [14, 12]. If we assume a
fiducial model with power spectra $C_{\ell}^{\mathrm{fid}}$, while the true
model that describes our Universe is $C_{\ell}^{\mathrm{true}}$, the bias
$\Delta\theta_{\alpha}$ to the parameter $\theta_{\alpha}$ is given by

$\Delta\theta_{\alpha}=\sum_{\beta}F^{-1}_{\alpha\beta}\sum_{\ell_{b},\ell_{b}^{\prime}}\frac{\partial
C_{\ell_{b}}^{\mathrm{fid}}}{\partial\theta_{\beta}}\mathbb{C}^{-1}_{\ell_{b}\ell_{b}^{\prime}}\left(C_{\ell_{b}^{\prime}}^{\mathrm{true}}-C_{\ell_{b}^{\prime}}^{\mathrm{fid}}\right).$
(20)

Here the Fisher matrix $F_{\alpha\beta}$ is calculated with the fiducial
model. We note that the parameter biases obtained using Eq. 20 are technically
conservative for the case where one has a prior on the value of a parameter
(see Appendix A). This is because, while Eq. 20 allows one to specify a
parameter uncertainty from prior data via the Fisher matrix $F_{\alpha\beta}$,
it does not allow one to specify the mean value of this prior, potentially
resulting in artificially larger bias estimates than is the actual case. The
bias estimates from Eq. 20 can be cross checked with an MCMC method to verify
consistency, as done in Appendix A.

###  IV.2 Likelihood and Markov Chain

We use an MCMC method to verify our Fisher parameter forecasts. With an MCMC
analysis, we want to recover the posterior probability distribution
$P(\vec{\theta}\ |\hat{d})$ for a set of model parameters $\vec{\theta}$ that
are most likely to describe the observed data $\hat{d}$. The posterior
distribution is proportional to the product of the likelihood function,
$\mathcal{L}(\hat{d}\ |\vec{\theta})$, and the prior probability distribution,
$P(\vec{\theta})$, for the parameter values $\vec{\theta}$. The posterior
distribution is estimated by drawing random samples of sets of parameters,
$\vec{\theta}$, via an MCMC sampler, and evaluating the likelihood and prior
probability for each [55].

For the results presented in this work, we use the mcmc sampler [56, 57] from
the Cobaya [17] package to sample from the likelihood functions defined in
Eqs. 16 and 17, and we use CAMB to calculate the theory
$C_{\ell}(\vec{\theta})$ or $f_{j}(\vec{\theta})$ at each step in parameter
space.212121Note that Cobaya calls CAMB itself, but we added an interface
within our likelihood to obtain the delensed spectra from CAMB via Cobaya. The
priors adopted for each MCMC analysis are listed in Table 2; all priors are
uniform except for the Gaussian prior on the reionization optical depth from
Planck of $\tau=0.054\pm 0.007$ [1].

MCMC runs in general converge faster if given an initial proposal matrix
consisting of the anticipated covariance between the cosmological parameters.
We generate such a proposal matrix by first using a Fisher matrix to obtain
the conditional posterior width for each parameter, given by
$1/\sqrt{F_{\alpha\alpha}}$ for parameter $\theta_{\alpha}$. We use this to
determine the step size of each parameter in the sampler. We then generate
mock data and run a quick preliminary MCMC using the default CAMB accuracy,
except for setting lens_potential_accuracy=1 to turn on the non-linear matter
power spectrum (see Appendix A).222222Increasing the CAMB accuracy settings
higher than this does not change the parameter error bars as long as the mock
data is generated at the same accuracy; however, we show in Appendix A that
low CAMB accuracy settings in general may result in a bias when running MCMC
chains on real data. We use the resulting MCMC chains to calculate a parameter
covariance matrix, which we use as a proposal matrix for subsequent MCMC runs
with our baseline CAMB accuracy.

Figure 7: Forecasted cosmological parameter constraints for a
$\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$ model from delensed CMB $TT$,
$TE$, $EE$, and $BB$ power spectra and the lensing power spectrum,
$\kappa\kappa$, for SO-like, CMB-S4-like, and CMB-HD-like surveys. The
experimental configurations for each survey are listed in Table 1. We also
include expected DESI BAO data [16] to these forecasts. The parameters that
see the most improvement from a CMB-HD type survey compared to precursor
experiments are $n_{\mathrm{s}}$ and $N_{\mathrm{eff}}$. We give the one-
dimensional marginalized parameter errors in Table 3.

We consider the MCMC parameter chains to be converged when they reach a value
of about $R-1\leq 0.01$ after removing the first 30% of each chain as “burn-
in”, where $R-1$ is the Gelman-Rubin convergence parameter [58]. We find that
the likelihood for delensed CMB-HD plus DESI BAO for a $\Lambda$CDM +
$N_{\mathrm{eff}}$ \+ $\sum m_{\nu}$ model takes about 12 hours on
NERSC232323https://www.nersc.gov/ using our baseline accuracy to reach
$R-1=0.03$, and about 16 hours to reach $R-1=0.015$. These times are obtained
after including a proposal matrix generated as described above. We confirm
that the marginalized mean parameter values that we obtain match the fiducial
values used to generate the mock data, listed in Table 2. We also confirm that
we recover parameter errors that match those obtained from the Fisher analysis
described above (see Appendix B and Fig. 18 for details). Since the MCMC run
takes about a day, after we verify a few base cases, we use the Fisher method
to explore variations of those cases (e.g. changes in $\ell_{\rm{max}}$,
inclusion of foregrounds or delensing, etc.).

##  V Parameter Forecasts

|  |  $1\sigma$ Error | Ratio of $1\sigma$ Errors
---|---|---|---
Parameter | Fiducial | SO + | CMB-S4 + | CMB-HD + | CMB-HD / SO | CMB-HD / CMB-S4
$\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$ |  | DESI BAO | DESI BAO | DESI BAO |  |
$\Omega_{\mathrm{b}}h^{2}$ |  $0.022370$ |  $0.000057$ |  $0.000039$ |  $0.000026$ |  $0.46$ | $0.67$
$\Omega_{\mathrm{c}}h^{2}$ |  $0.12000$ |  $0.00074$ |  $0.00056$ |  $0.00041$ |  $0.55$ | $0.73$
$\ln(10^{10}A_{\mathrm{s}})$ |  $3.044$ |  $0.012$ |  $0.011$ |  $0.0098$ |  $0.84$ | $0.86$
$n_{\mathrm{s}}$ |  $0.9649$ |  $0.0030$ |  $0.0025$ |  $0.0013$ |  $0.43$ | $0.52$
$\tau$ |  $0.0544$ |  $0.0061$ |  $0.0060$ |  $0.0052$ |  $0.85$ | $0.87$
$H_{0}$ [km s-1 Mpc-1] |  $67.36$ |  $0.36$ |  $0.32$ |  $0.29$ |  $0.81$ | $0.91$
$N_{\mathrm{eff}}$ |  $3.046$ |  $0.043$ |  $0.030$ |  $0.014$ |  $0.33$ | $0.47$
$\sum m_{\nu}$ [eV] |  $0.06$ |  $0.030$ |  $0.029$ |  $0.025$ |  $0.83$ | $0.86$
Table 3: Forecasted cosmological parameter constraints for a $\Lambda$CDM + $N_{\rm{eff}}+\Sigma m_{\nu}$ model from SO-like, CMB-S4-like, and CMB-HD-like delensed $TT$, $TE$, $EE$, $BB$ and $\kappa\kappa$ power spectra when combined with DESI BAO. The first two columns list the parameters and their fiducial values. The following three columns list the forecasted marginalized $1\sigma$ uncertainties on each parameter for the three experiments, and the last two columns list the ratios of these values. The forecasts for CMB-HD include expected residual extragalactic foregrounds in the temperature power spectrum. All forecasts include a $\tau$ prior of $\tau=0.054\pm 0.007$ from Planck [1]. We find improvement in all parameters considered for a CMB-HD-like experiment compared to precursor surveys; we also find the most significant improvement for $n_{\mathrm{s}}$ and $N_{\mathrm{eff}}$, of about a factor of two or more. These forecasts are also depicted in Fig. 7.  |  | $\ell_{\mathrm{max}},\leavevmode\nobreak\ L_{\mathrm{max}}=20,000$ | $\ell_{\mathrm{max}},L_{\mathrm{max}}=10,000$
---|---|---|---
Parameter | Fiducial | HD Lensed | HD Lensed | HD Delensed |  HD Delensed | HD Delensed
$\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$ |  |  | \+ FG | \+ FG |  \+ FG + DESI BAO | \+ FG + DESI BAO
$\Omega_{\mathrm{b}}h^{2}$ |  $0.022370$ |  $0.000032$ |  $0.000033$ |  $0.000027$ |  $0.000026$ | $0.000026$
$\Omega_{\mathrm{c}}h^{2}$ |  $0.12000$ |  $0.00064$ |  $0.00065$ |  $0.00058$ |  $0.00041$ | $0.00041$
$\ln(10^{10}A_{\mathrm{s}})$ |  $3.044$ |  $0.011$ |  $0.011$ |  $0.011$ |  $0.0098$ | $0.010$
$n_{\mathrm{s}}$ |  $0.9649$ |  $0.0023$ |  $0.0024$ |  $0.0021$ |  $0.0013$ | $0.0014$
$\tau$ |  $0.0544$ |  $0.0056$ |  $0.0057$ |  $0.0054$ |  $0.0052$ | $0.0054$
$H_{0}$ [km s-1 Mpc-1] |  $67.36$ |  $0.72$ |  $0.74$ |  $0.65$ |  $0.29$ | $0.29$
$N_{\mathrm{eff}}$ |  $3.046$ |  $0.017$ |  $0.018$ |  $0.015$ |  $0.014$ | $0.015$
$\sum m_{\nu}$ [eV] |  $0.06$ |  $0.050$ |  $0.051$ |  $0.047$ |  $0.025$ | $0.026$
Table 4: Forecasted cosmological parameter constraints from CMB-HD $TT$, $TE$,
$EE$, $BB$ and $\kappa\kappa$ power spectra for a $\Lambda$CDM +
$N_{\rm{eff}}+\Sigma m_{\nu}$ model. All forecasts include a $\tau$ prior of
$\tau=0.054\pm 0.007$ from Planck [1]. The first two columns list the
parameters and their fiducial values. The following two columns list their
forecasted marginalized $1\sigma$ uncertainties when using lensed spectra with
or without foregrounds in the temperature maps. The fifth column shows the
forecast when delensing the CMB spectra, and the sixth column shows the change
when including DESI BAO data [16]. In each of these cases we use CMB and CMB
lensing multipoles out to $\ell_{\mathrm{max}},L_{\mathrm{max}}=20,000$ for
CMB-HD. The last column lists the same information as the sixth column, but
instead using a maximum multipole of
$\ell_{\mathrm{max}},L_{\mathrm{max}}=10,000$ for both CMB and CMB lensing
power spectra. We see in particular that $n_{\mathrm{s}}$ and $N_{\rm{eff}}$
constraints are tightened when including multipoles beyond 10,000 for a CMB-HD
type survey. We show corresponding forecasts for $\Lambda$CDM and $\Lambda$CDM
+ $N_{\rm{eff}}$ models in Appendix C.

For the discussion below, we show parameter constraints using Fisher matrix
estimation since it provides a faster way to explore the impact of several
effects. However, as mentioned above, we confirm that our parameter forecasts
are consistent when using either a likelihood and MCMC or Fisher matrix
estimation for a subset of cases (see Appendix B for details). For the SO-like
and S4-like experiments, we also confirm that our parameter constraints are
consistent with those forecasted in [4] and [5], respectively.

In Fig. 7, we show the cosmological parameter constraints for the
$\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$ model from delensed CMB $TT$,
$TE$, $EE$, $BB$ spectra and lensing $\kappa\kappa$ spectra when combined with
DESI BAO data. For the CMB-HD-like experiment (red), we have included residual
extragalactic foregrounds in the temperature data. For comparison, we forecast
parameter constraints for S4-like (green) and SO-like (yellow) experiments. We
show the $1\sigma$ and $2\sigma$ parameter contours as the dark and lighter
colors, respectively, and we list the $1\sigma$ marginalized parameter
constraints for all three experiments in Table 3. We find improvement in all
parameters considered for a CMB-HD-like experiment compared to precursor
surveys; we also find the most significant improvement for $n_{\mathrm{s}}$
and $N_{\mathrm{eff}}$, of about a factor of two or more (see last two columns
of Table 3). In particular, for CMB-HD, we find
$\sigma(n_{\mathrm{s}})=0.0013$ and $\sigma(N_{\mathrm{eff}})=0.014$.

In the following sections, we focus on a CMB-HD-like experiment and discuss
the effects of residual extragalactic foregrounds, delensing of the acoustic
peaks, and the combination with DESI BAO data on cosmological parameter
constraints. The resulting parameter uncertainties in each case for a
$\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$ model are listed in Table 4; we
present results for a $\Lambda$CDM and a $\Lambda$CDM+$N_{\mathrm{eff}}$ model
in Appendix C. In the last column of Tables 4, 7, and 8, we examine the effect
of using a lower maximum multipole of
$\ell_{\mathrm{max}},L_{\mathrm{max}}=10,000$ and find, in particular, that
$n_{\mathrm{s}}$ and $N_{\rm{eff}}$ constraints are tightened when including
multipoles beyond 10,000.

Figure 8: Shown is the impact of including expected residual foregrounds (FG)
on the forecasted constraints for a $\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum
m_{\nu}$ model from CMB-HD lensed $TT$, $TE$, $EE$, $BB$ and $\kappa\kappa$
power spectra. Here we show a subset of parameter constraints for the Hubble
constant $H_{0}$, the effective number of relativistic species
$N_{\mathrm{eff}}$, and the scalar spectral index $n_{\mathrm{s}}$. We see
that the expected foregrounds do not significantly increase parameter
uncertainties; we emphasize that while these are expected foreground levels
from preliminary analyses [15], it remains to be demonstrated with realistic
simulations that such foreground levels can be achieved in practice.

###  V.1 Impact of Residual Foregrounds

We examine the effect of including residual extragalactic foregrounds in
temperature maps, as described in Section II.3, on the parameter constraints
for CMB-HD from lensed CMB $TT$, $TE$, $EE$, and $\kappa\kappa$ power spectra.
In Fig. 8, we show the parameter constraints without (orange) and with (green)
foregrounds for a subset of parameters from a
$\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$ model; the full set of
forecasted parameter errors are listed in the third and fourth columns of
Table 4. We see that the residual foregrounds do not significantly increase
the errors on any of the parameters considered here. The ultra-high resolution
and low noise of CMB-HD is critical to reducing the foreground levels to those
shown in the top panel of Fig. 1.

###  V.2 Impact of Delensing

Figure 9: Shown are forecasted constraints for a
$\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$ model from lensed (green) and
delensed (blue) CMB-HD $TT$, $TE$, $EE$, $BB$ and $\kappa\kappa$ power
spectra. Here we show a subset of parameter constraints for $H_{0}$,
$N_{\mathrm{eff}}$, and $n_{\mathrm{s}}$, and show the full parameter
constraints in Fig. 11. We see that delensing the CMB spectra reduces the
parameter uncertainties.

In Fig. 9, we show the parameter constrains from CMB-HD lensed (green) and
delensed (blue) CMB spectra and the lensing spectrum for a subset of
parameters from a $\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$ model. In
both cases, we include residual extragalactic foregrounds. We see smaller
uncertainties in both $n_{\mathrm{s}}$ and $N_{\mathrm{eff}}$ after delensing.
We show the comparison between lensed and delensed forecasts for the full
parameter set in Fig. 11. The fifth column of Table 4 gives the marginalized
$1\sigma$ parameter errors when including delensing, and shows considerable
improvement for all the parameters considered, even when also including
residual foregrounds. We attribute this to the removal of much of the off-
diagonal covariance on scales where CMB-HD can delens efficiently, as shown in
the right panel of Fig. 6. Tables 7 and 8 show the impact of delensing for
$\Lambda$CDM and $\Lambda$CDM+$N_{\mathrm{eff}}$ models, respectively. In
Table 9, we show that even when including DESI BAO, delensing still improves
parameter errors, especially for $N_{\mathrm{eff}}$.

###  V.3 Impact of DESI BAO

Figure 10: Shown are the forecasted constraints for a
$\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$ model from CMB-HD delensed
$TT$, $TE$, $EE$, $BB$ and $\kappa\kappa$ power spectra, with and without
including expected DESI BAO data [16] (blue and red, respectively). Here we
show a subset of parameter constraints for $H_{0}$, $N_{\mathrm{eff}}$, and
$n_{\mathrm{s}}$, and show the full parameter constraints in Fig. 11. We see
that the inclusion of DESI BAO data reduces uncertainties on $H_{0}$ and
$n_{\mathrm{s}}$, but has less impact on $N_{\mathrm{eff}}$.

We show in Figs. 10 and 11 the parameter constraints for CMB-HD delensed
$TT,TE,EE,BB$ and $\kappa\kappa$ spectra, without (blue) and with (red) DESI
BAO data. We see that adding BAO data tightens all the parameter
uncertainties, and most significantly impacts $\Omega_{\mathrm{c}}h^{2}$,
$H_{0}$, $n_{\mathrm{s}}$, and $\sum m_{\nu}$. The sixth column of Tables 4,7
and 8 give the marginalized $1\sigma$ parameter errors when including DESI BAO
for $\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$, $\Lambda$CDM, and
$\Lambda$CDM+$N_{\mathrm{eff}}$ models, respectively.

Figure 11: Shown are the forecasted constraints for a
$\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$ model for CMB-HD. Here we show
constraints from lensed CMB $TT$, $TE$, $EE$, $BB$ and $\kappa\kappa$ power
spectra (green), delensed CMB $TT$, $TE$, $EE$, $BB$ and $\kappa\kappa$ power
spectra (blue), and delensed CMB $TT$, $TE$, $EE$, $BB$ and $\kappa\kappa$
power spectra plus DESI BAO [16] (red). We see that delensing the CMB spectra
and adding DESI BAO data reduce uncertainties on cosmological parameters.

###  V.4 Impact of Baryonic Physics

The low noise and high resolution of CMB-HD allows it to be sensitive to small
scales ($k>0.5\leavevmode\nobreak\ h$Mpc-1) where the effects of baryonic
feedback on the matter distribution become important. On these scales,
feedback from active galactic nuclei (AGN) emissions, for example, can push
mass out from the centers of halos, while baryonic cores, formed from star-
formation and cooling, concentrate mass in the center of halos. These effects
shift the mass distribution, which can be measured with CMB lensing.

In addition, the hot gas traced by the thermal and kinetic Sunyaev-Zel’dovich
effects (kSZ and tSZ, respectively) is also redistributed by baryonic
feedback. This makes the kSZ and tSZ effects important external measurements
of the amount of feedback that has occurred and the impact of that feedback.
In particular, cross-correlations between CMB lensing and the kSZ/tSZ effects
can provide powerful additional constraints on the behavior of baryonic
effects [59].

To model the impact of baryonic feedback, we use the updated HMcode-2020 model
from [60], and in particular adopt their single-parameter baryonic feedback
model characterized by the parameter
$\mathrm{log}_{10}({T_{\mathrm{AGN}}}/{\mathrm{K}})$242424Previous versions of
HMCode did not include star-formation in their baryonic model, and only
predicted power suppression [60]. The HMcode-2020 model has six free
parameters, but [60] find a single-parameter variant is only slightly less
accurate (e.g. see Fig 5 of [60]).; we use a fiducial value of
$\mathrm{log}_{10}({T_{\mathrm{AGN}}}/{\mathrm{K}})=7.8$ and a step size of
$\pm 0.05$ in the Fisher analysis. This model was used in a recent analysis
combining KiDS-1000 optical lensing data with tSZ measurements from Planck and
ACT. This analysis measured the cross-correlation between the lensing shear
and tSZ maps, and measured
$\mathrm{log}_{10}({T_{\mathrm{AGN}}}/{\mathrm{K}})$ with a $1\sigma$
uncertainty better than $6\%$ [59].

We show in Fig. 12 the impact of freeing this baryonic feedback parameter on a
subset of parameters for a $\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum
m_{\nu}$+baryonic feedback model from CMB-HD delensed $TT$, $TE$, $EE$, $BB$
and $\kappa\kappa$ power spectra plus DESI BAO (light blue). We also show the
marginalized $1\sigma$ parameter errors in the second column of Table 5. We
see that the parameter errors do not increase substantially when freeing this
baryonic parameter. We also see that CMB-HD can constrain this baryonic
feedback parameter to an accuracy of 0.45% without including any information
from kSZ or tSZ. This tight constraint on baryonic feedback is likely due to
the precision CMB-HD will have in measuring the small-scale lensing power
spectrum.

We also add a prior on the baryonic feedback parameter
$\mathrm{log}_{10}({T_{\mathrm{AGN}}}/{\mathrm{K}})$ that we can expect from
CMB-HD measurements of the tSZ and kSZ effects in cross-correlation with CMB
lensing. We estimate this prior by extrapolating from the 6% constraint
measured by [59] using KiDS-1000 lensing maps cross-correlated with Planck
plus ACT tSZ maps over 1000 square degrees of sky. We note that CMB-HD will
survey about 24,000 square degrees of sky, gaining about a factor of 5
improvement over the current feedback constraint. In addition, the CMB-HD
temperature map noise level will be over 20 times deeper than that of Planck
plus ACT. Moreover, the CMB-HD lensing map will be a factor of a few higher in
signal-to-noise ratio than that of KiDS-1000, and we can anticipate also
folding in information from the cross-correlation of kSZ maps with CMB lensing
(not included in the KiDS-1000 analysis). To be conservative, we assume an
overall improvement of two orders of magnitude in ability to constrain
$\mathrm{log}_{10}({T_{\mathrm{AGN}}}/{\mathrm{K}})$ from CMB-HD tSZ and kSZ
cross-correlations with CMB lensing compared to the current constraint. (As
mentioned above, CMB-HD already can improve the constraint on
$\mathrm{log}_{10}({T_{\mathrm{AGN}}}/{\mathrm{K}})$ by an order of magnitude
without any SZ data.) Thus we apply a 0.06% prior on this baryonic feedback
parameter, and show the results in Fig. 12 (dashed black) and Table 5 (third
column). We see that the inclusion of this additional prior from SZ data
returns the parameter errors to what they were before freeing the feedback
parameter (shown in the first column of Table 5).

Figure 12: Shown is the impact of freeing the baryonic feedback parameter,
$\mathrm{log}_{10}({T_{\mathrm{AGN}}}/{\mathrm{K}})$, discussed in Section
V.4, on the forecasted constraints for a $\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum
m_{\nu}$+baryonic feedback model (cyan) from CMB-HD delensed $TT$, $TE$, $EE$,
$BB$ and $\kappa\kappa$ power spectra plus DESI BAO. We also show the
parameter constraints when fixing this parameter (red) or adding a prior on
the baryonic feedback parameter (dashed black); we anticipate such a prior can
be obtained from a joint analysis of CMB-HD kSZ, tSZ, and lensing data (see
Section V.4 for details). Here we show a subset of parameter constraints for
$H_{0}$, $N_{\mathrm{eff}}$, and $n_{\mathrm{s}}$ from a Fisher analysis.

We use Eq. 20 in Section IV.1 to predict the bias on the estimated
cosmological parameters when AGN feedback is neglected in the model. In this
case, we allow $C_{\ell}^{\mathrm{fid}}$ to be the CDM-only model whereas the
true model $C_{\ell}^{\mathrm{true}}$ is CDM+AGN feedback with
$\mathrm{log}_{10}({T_{\mathrm{AGN}}}/{\mathrm{K}})=7.8$. We also use the
combination of CMB-HD and DESI BAO data to calculate the Fisher matrix,
$F_{\alpha\beta}$, in Eq. 20. In Fig. 13, we show the forecasted parameter
biases from the combination of lensed (purple) or delensed (blue) CMB-HD and
DESI BAO data when assuming a fiducial CDM-only model, as opposed to the true
CDM+feedback model. We can see in Fig. 13 that, while delensing reduces the
bias due to an incorrect baryonic feedback model (as discussed in [14, 61,
12]), for many parameters the remaining bias is much larger than the
statistical $1\sigma$ error.

We also show in Fig. 13 the result of an MCMC run when marginalizing over the
baryonic feedback parameter
$\log_{10}\left(T_{\mathrm{AGN}}/\mathrm{K}\right)$ discussed above, applying
only a uniform prior of
$\log_{10}\left(T_{\mathrm{AGN}}/\mathrm{K}\right)\in[7.6,\leavevmode\nobreak\
8.0]$ as suggested by [60]. We see that this marginalization over feedback
models removes the parameter biases at the expense of some constraining power.
However, as we show in Fig. 12, adding a prior on the baryonic feedback
parameter anticipated from CMB-HD SZ measurements can mitigate this increase
in parameter error. We note that the MCMC run with the baryonic feedback
parameter free takes about the same time to converge on NERSC as a run with
the feedback parameter fixed.

Parameter | CDM | \+ feedback | \+ SZ prior
---|---|---|---
$\Omega_{\mathrm{b}}h^{2}$ |  $0.000026$ |  $0.000028$ | $0.000027$
$\Omega_{\mathrm{c}}h^{2}$ |  $0.00041$ |  $0.00046$ | $0.00041$
$\ln(10^{10}A_{\mathrm{s}})$ |  $0.0095$ |  $0.010$ | $0.010$
$n_{\mathrm{s}}$ |  $0.0013$ |  $0.0021$ | $0.0013$
$\tau$ |  $0.0051$ |  $0.0056$ | $0.0056$
$100\theta_{\mathrm{MC}}$ |  $0.000058$ |  $0.000064$ | $0.000060$
$N_{\mathrm{eff}}$ |  $0.014$ |  $0.020$ | $0.014$
$\sum m_{\nu}$ [eV] |  $0.024$ |  $0.026$ | $0.027$
$\log_{10}\left(T_{\mathrm{AGN}}/\mathrm{K}\right)$ | — |  $0.040$ | $0.0047$
$H_{0}$ [km s-1 Mpc-1] |  $0.28$ |  $0.29$ | $0.27$
$\sigma_{8}$ |  $0.0033$ |  $0.0033$ | $0.0033$
Table 5: Shown are the cosmological constraints from the combination of CMB-HD
delensed $TT$, $TE$, $EE$, $BB$, and $\kappa\kappa$ power spectra with DESI
BAO for the baseline $\Lambda$CDM + $N_{\mathrm{eff}}$ \+ $\sum m_{\nu}$
model, a $\Lambda$CDM + $N_{\mathrm{eff}}$ \+ $\sum m_{\nu}$ \+ baryonic
feedback model from HMCode-2020 [60], and a $\Lambda$CDM + $N_{\mathrm{eff}}$
\+ $\sum m_{\nu}$ \+ baryonic feedback model including a 0.06% prior on the
feedback parameter $\mathrm{log}_{10}({T_{\mathrm{AGN}}}/{\mathrm{K}})$
expected from from a joint analysis of CMB-HD kSZ, tSZ, and lensing data (see
Section V.4 for details). All results shown here are from the likelihood and
MCMC chains as opposed to Fisher forecasts. We also include
$100\theta_{\mathrm{MC}}$, and separate the two derived parameters, $H_{0}$
and $\sigma_{8}$, which is the linear root-mean-square matter fluctuations
today. Figure 13: Here we show the expected bias to the parameter constraints
from a combination of lensed or delensed (purple and blue, respectively) CMB-
HD $TT$, $TE$, $EE$, $BB$ and $\kappa\kappa$ power spectra plus DESI BAO, when
baryonic feedback effects are neglected. These biases were computed from Eq.
20, using a CDM-only model as the fiducial model, and assuming the true model
includes the effects of baryonic feedback. We center the CDM-only contours at
the estimated biased parameter values, and take the error ellipses from the
Fisher estimates. In cyan we show the parameter constraints for delensed CMB-
HD plus DESI BAO data from an MCMC run where we marginalize over the feedback
parameter $\log_{10}\left(T_{\mathrm{AGN}}/\mathrm{K}\right)$ discussed in
Section V.4. We see that marginalization over feedback models removes the
parameter biases.

##  VI Discussion

In this work, we present the parameter forecasts for a CMB-HD survey. We
contrast these forecasts with precursor CMB experiments, showing that the
lower noise and higher multipoles of a CMB-HD survey can lead of significant
improvement in many of the cosmological parameters, most notably in the scalar
spectral index, $n_{\mathrm{s}}$, and in the effective number of light
relativistic species, $N_{\mathrm{eff}}$. Specifically we find for delensed
CMB-HD $TT,TE,EE,BB$ and $\kappa\kappa$ spectra plus DESI BAO:

$\sigma\left(n_{\mathrm{s}}\right)=0.0013\quad(68\%,\leavevmode\nobreak\
TT,TE,EE,BB,\kappa\kappa+\mathrm{BAO}),$
$\sigma\left(N_{\mathrm{eff}}\right)=0.014\quad(68\%,\leavevmode\nobreak\
TT,TE,EE,BB,\kappa\kappa+\mathrm{BAO}).$

We find that delensing all the CMB spectra as well as including DESI BAO data
are both necessary to achieve these tight constraints. We also find that
including multipoles out to $\ell_{\mathrm{max}},L_{\mathrm{max}}=20,000$ is
also required. Since these multipoles are well into the non-linear regime, we
must include baryonic effects in modelling the lensing in both the CMB power
spectra and CMB lensing spectra. We find that marginalizing over baryonic
effects can mitigate potential bias in parameters at the expense of some
constraining power. However, tSZ and kSZ measurements by CMB-HD offer an
independent handle on baryonic feedback effects, and folding in that
additional information can effectively eliminate the increase in parameter
errors due to uncertainty in the baryonic physics.

The $N_{\mathrm{eff}}$ parameter uncertainty achieved by CMB-HD+DESI is
particularly interesting since any new light particle that was in thermal
equilibrium at any time after the Universe reheated must change
$N_{\mathrm{eff}}$ by at least $\Delta N_{\mathrm{eff}}\geq 0.027$ [62, 63].
Here reheating refers to the end of inflation and the beginning of the “Big
Bang”. With $\sigma\left(N_{\mathrm{eff}}\right)=0.014$, CMB-HD+DESI can
either rule out or detect any new light particle species with at least 95%
confidence.

As a specific example of why this $N_{\mathrm{eff}}$ constraint is valuable,
we consider the QCD axion, which is a well-motivated candidate for being the
dark matter and solving the Strong CP problem [64, 65, 66]. We note that if
the reheating temperature of the Universe is high enough that the QCD axion
thermalized, the $N_{\mathrm{eff}}$ constraint above can potentially rule out
the QCD axion in a model-independent way, or lead to a detection.

We re-plot Fig. 3 of [67] in Fig. 14 to highlight the QCD axion masses,
$m_{\phi}$, that would be ruled out as a function of reheating temperature,
$T_{\mathrm{R}}$, if CMB-HD+DESI finds no increase in $N_{\mathrm{eff}}$ at
the 95% confidence level. To forecast these constraints, we use Eq. 2.11 of
[67] to calculate the upper limit of the QCD axion coupling $g_{d}$ for a
given reheating temperature,

$g_{d}<1.3\times 10^{-14}\leavevmode\nobreak\
\mathrm{GeV}^{-2}\left(\frac{T_{\mathrm{R}}}{10^{10}\leavevmode\nobreak\
\mathrm{GeV}}\right)^{-1/2},$ (21)

and use Eq. 16 of [68] to relate this to the coupling constant $f_{a}$ via

$g_{d}\approx\frac{2.4\times 10^{-16}\leavevmode\nobreak\
\mathrm{e}\leavevmode\nobreak\ \mathrm{cm}}{f_{a}},$ (22)

where $1\leavevmode\nobreak\ \mathrm{cm}^{-1}=1.97\times
10^{-14}\leavevmode\nobreak\ \mathrm{GeV}$ and $\mathrm{e}=0.3$. We then
relate the coupling constant to the axion mass using

$m_{\phi}=0.6\leavevmode\nobreak\
\mathrm{eV}\left(\frac{10^{7}\leavevmode\nobreak\ \mathrm{GeV}}{f_{a}}\right)$
(23)

from [64, 69, 70]. Given the well-motivated nature of the QCD axion and the
many efforts to detect it underway, it is worth highlighting this model-
independent additional approach to probing this new physics.

We also note that the $n_{\mathrm{s}}$ constraint above is about a factor of
two tighter than from precursor CMB surveys. While considerable attention has
been devoted to improving constraints on primordial gravitational waves via
the tensor-to-scalar ratio, $r$, constraining the scalar spectral index can
also rule out interesting inflationary scenarios [5].

Figure 14: Here we show constraints on the QCD axion mass from delensed CMB-HD
plus DESI BAO data. The shaded region shows the region of axion mass,
$m_{\phi}$, and reheating temperature, $T_{\mathrm{R}}$, that would be
excluded if any measurement rules out $\Delta N_{\mathrm{eff}}\geq 0.027$,
which we show in this work CMB-HD+DESI BAO does at the 95% confidence level.

##  VII Conclusion

While we present parameter forecasts for CMB-HD above, we note that there are
several ways that these forecasts can be improved and made more robust. As we
have stressed a few times, the most challenging aspect of a CMB-HD survey will
be removing and mitigating the impact of extragalactic foregrounds.
Demonstrating that this can be achieved with realistic end-to-end simulations
is an area of ongoing research, and will be critical for achieving the science
presented here. In addition, more optimal lensing estimators have been
proposed and are being explored [36, 71, 37, 9, 33, 72, 73, 38, 74], with the
potential to yield higher lensing signal-to-noise ratios than assumed in this
work, as well as to provide better paths towards foreground immunity. We also
leave to future work exploring the potential gain in lensing signal-to-noise
ratio from exploiting the higher-order lensing corrections contained in the N1
signal. Furthermore, we have only explored cold dark matter models in this
work, with and without baryonic feedback, but note that CMB-HD has sensitivity
to alternate dark matter models as well, which will be investigated in
subsequent work.

We make the Fisher estimation code used here public along with a Jupyter
notebook detailing how to generate new Fisher derivatives for different
models. We also make the likelihood code public and integrate it with Cobaya
using CAMB. We hope that this will facilitate additional cosmological
parameter forecasts for a CMB-HD survey.

###### Acknowledgements.

The authors thank Dongwon Han, Gil Holder, Selim Hotinli, Mathew
Madhavacheril, Joel Meyers, Vivian Miranda, Rugved Pund, Cynthia Trendafilova,
and Alexander van Engelen for useful discussions. NS also thanks Itay Bloch,
Rouven Essig, Peter Graham, Maxim Pospelov, Mauro Valli, and the SCGP Lighting
New Lampposts Workshop for useful discussions about the QCD axion. AM, NS, and
MR acknowledge support from DOE award number DE-SC0020441. AM and NS also
acknowledge support from NSF award number 1907657. This research used
resources of the National Energy Research Scientific Computing Center (NERSC),
a U.S. Department of Energy Office of Science User Facility located at
Lawrence Berkeley National Laboratory, operated under Contract No. DE-
AC02-05CH11231 using NERSC award HEP-ERCAPmp107. NS also acknowledges the
Aspen Center for Physics, which is supported by National Science Foundation
grant PHY-1607611, for supporting a workshop on new discoveries in the era of
high-resolution, low-noise CMB experiments, which stimulated useful
discussions.

##  Appendix A Einstein-Boltzmann Code Accuracy Requirements

As discussed in [14], future high-resolution and low noise CMB experiments
will require theory calculations for power spectra beyond the standard default
CAMB accuracy settings for the Einstein-Boltzmann solver. To determine the
accuracy settings required for CMB-HD, we focus on the parameter bias induced
by modelling the expected data with accuracy settings that are lower than
necessary. For this we use Eq. 20 discussed in Section IV.1, and we explore
four accuracy parameters provided by CAMB (i.e. lens_potential_accuracy,
AccuracyBoost, lAccuracyBoost, and lSampleBoost). We also by default set
NonLinear=model.NonLinear_both to enable the non-linear calculation of both
the matter power spectrum and the lensing spectrum. We note that as long as
non-linear corrections are included, changing the four CAMB accuracy settings
above does not change the error bars obtained from either Fisher or MCMC
analyses, as long as the mock data is generated with the same accuracy.

Figure 15: Here we compare the lensing power spectrum for different values of
the CAMB accuracy parameter lens_potential_accuracy. The top panel shows the
fractional difference between $C_{L}^{\kappa\kappa}$ and
$C_{L}^{\kappa\kappa,\mathrm{true}}$, while the bottom panel shows the
difference as a fraction of the CMB-HD error bar on the lensing power at each
multipole, $\sigma_{L}^{\mathrm{HD}}$. $C_{L}^{\kappa\kappa,\mathrm{true}}$
uses the accuracy settings we describe in the text for our approximation of
the “true Universe”. We see that the lensing power spectrum converges to the
“true Universe” model when lens_potential_accuracy is increased, and for
lens_potential_accuracy > 30 the difference between the spectra is less than
10% of the CMB-HD error bar for a given $L$. Here we fix the other parameters
to AccuracyBoost = 1.1, lAccuracyBoost = 3.0, and lSampleBoost = 3.0. Since
the “true Universe” model has AccuracyBoost = 3.0, lAccuracyBoost = 5.0, and
lSampleBoost = 5.0, we see that varying these other parameters within these
ranges has minimal impact given the CMB-HD error bars.

We find that the lens_potential_accuracy setting has the largest impact on
parameter biases. We show in the top panel of Fig. 15, the change in the
lensing power spectrum for different values of lens_potential_accuracy, noting
that for higher accuracy settings most of the change is at high lensing
multipoles. We show in the bottom panel of Fig. 15, the change in the lensing
power spectrum divided by the expected CMB-HD uncertainty per lensing
multipole. We find that above lens_potential_accuracy = 30, the lensing power
spectrum theory curves converge with a difference smaller than a tenth of the
CMB-HD error bar per multipole. Thus we take lens_potential_accuracy = 40 to
be the “true Universe” model for this work, and compare our results to this
model.

We find that after setting lens_potential_accuracy = 40, varying the
lAccuracyBoost and lSampleBoost parameters results in parameter biases that
are well below $1\sigma$ of the expected CMB-HD+BAO parameter errors for a
$\Lambda$CDM + $N_{\mathrm{eff}}$ \+ $\sum m_{\nu}$ model. The default CAMB
settings for these parameters are 1.0, and we increase the settings to 5.0 for
each for our “true Universe” model to be conservative.

Figure 16: Here we show the expected bias on each cosmological parameter as a
fraction of its $1\sigma$ error for CMB-HD+DESI when we calculate the delensed
CMB and lensing power spectra with different values of the CAMB accuracy
parameters lens_potential_accuracy and AccuracyBoost. The other parameters are
fixed to values of lAccuracyBoost = 3.0 and lSampleBoost = 3.0, with the
“true” model in Eq. 20 calculated with the accuracy settings AccuracyBoost =
3.0, lAccuracyBoost = 5.0, and lSampleBoost = 5.0 (see text). We see that the
bias converges to below $0.3\sigma$ at lens_potential_accuracy = 30, and that
increasing the AccuracyBoost does not significantly decrease the bias.

In addition, we find that varying the AccuracyBoost parameter, also minimally
impacts parameter biases, and also results in biases well below $1\sigma$ of
the expected CMB-HD+BAO parameter errors. The CAMB default setting is
AccuracyBoost = 1.0, and we increase this to AccuracyBoost = 3.0 for our “true
Universe” model.

We show in Fig. 16, the bias on each parameter divided by the expected CMB-
HD+BAO $1\sigma$ parameter error for a $\Lambda$CDM + $N_{\mathrm{eff}}$ \+
$\sum m_{\nu}$ model. We use Eq. 20 to calculate each parameter bias, using
the “true Universe” model described above as $C_{\ell}^{\mathrm{true}}$. We
see that lens_potential_accuracy = 30 yields parameter biases that are less
than 0.3$\sigma$. We also find that varying the AccuracyBoost has minimal
impact.

In addition, we find that varying the AccuracyBoost has the most significant
impact on the computation time. We show in Fig. 17 the rapid increase in
computation time with increasing AccuracyBoost. Since we see from Fig. 16 that
an AccuracyBoost = 1.1 yields similar biases to a value double that, we choose
as our baseline setting AccuracyBoost = 1.1; we note that AccuracyBoost = 1.0
is the CAMB default, and we find a significant reduction in bias for a small
increase above that. Similarly, we find little reduction in bias for values of
lAccuracyBoost and lSampleBoost above 3.0.

Figure 17: Here we show the increase in the total CAMB computation time to
compute the lensed and delensed CMB spectra plus the CMB lensing power
spectrum for a maximum $\ell/L$ of 20,000 when the AccuracyBoost parameter is
increased. Here the other accuracy parameters are fixed to lAccuracyBoost =
3.0, lSampleBoost = 3.0, and lens_potential_accuracy = 30. The computation
times were measured on a NERSC Perlmutter login node. The red star indicates
the baseline accuracy settings used throughout this work, which shows that all
the spectra mentioned above can be calculated in a total time of 11 seconds.

Thus, we use as our baseline CAMB settings for CMB-HD in this work:

⬇

import camb

import numpy as np

lmax = 20100

pars = camb.CAMBparams()

pars.set_cosmology(H0=67.36, ombh2=0.02237, omch2=0.1200, tau=0.0544,
num_massive_neutrinos=1, mnu=0.06, nnu=3.046)

pars.InitPower.set_params(As=np.exp(3.044) * 1e-10, ns=0.9649)

pars.set_for_lmax(int(lmax)+500, lens_potential_accuracy=30, lens_margin=2050)

pars.set_accuracy(AccuracyBoost=1.1, lSampleBoost=3.0, lAccuracyBoost=3.0,
DoLateRadTruncation=False)

pars.NonLinear = camb.model.NonLinear_both

pars.NonLinearModel.set_params("mead2016")

We find that the computation of the lensed and delensed $TT,TE,EE,BB$ CMB
power spectra plus the CMB lensing power spectrum takes 11 seconds total to
run on the NERSC Perlmutter machine using the CAMB accuracy settings above.

We also run an MCMC to verify the accuracy of the Fisher estimated biases for
our baseline settings. We do this by generating mock data at the accuracy
settings of the “true Universe” model described above, and running the MCMC
with our baseline accuracy settings. We find consistency between the two
methods of bias estimation, and confirm that the parameter biases are well
below $1\sigma$ of the expected CMB-HD+BAO parameter errors.

##  Appendix B Consistency of Fisher and MCMC Methods

To verify that our Fisher forecasts accurately predict the expected parameter
errors, we create a likelihood and run MCMC chains. In Fig. 18 we show as
solid lines/shaded contours our MCMC run for the baseline CMB-HD delensed
$TT$, $TE$, $EE$, $BB$ and $\kappa\kappa$ power spectra (including
foregrounds), combined with mock DESI BAO data. We see that the MCMC recovers
the input values of the data well, which are indicated by the grey dashed
lines. We also overlay the results of the Fisher estimation method as red
dotted contours. We find good agreement between the two methods. We also
compare the $1\sigma$ parameter errors from each method in Table 6, finding
good consistency. Given the consistency of both methods, we use the Fisher
method for all parameter forecasts in this work, unless stated otherwise.

Figure 18: Here we show the parameter constraints for a $\Lambda$CDM +
$N_{\mathrm{eff}}$ \+ $\sum m_{\nu}$ model obtained from delensed CMB-HD $TT$,
$TE$, $EE$, $BB$ and $\kappa\kappa$ power spectra plus DESI BAO. The Fisher
forecasts are shown in red (dotted lines), and the results from running MCMC
chains are shown in blue (solid lines/shaded contours). The parameter values
used to generate the data are indicated by the grey dashed lines. We find good
agreement between parameter error estimates based on Fisher and MCMC analyses.

##  Appendix C Additional Parameter Forecasts

Below we show parameter constraints for CMB-HD in a $\Lambda$CDM model with
fixed $N_{\mathrm{eff}}=3.046$ and $\sum m_{\nu}=0.06$ eV (Table 7), and in a
$\Lambda$CDM + $N_{\mathrm{eff}}$ model with fixed $\sum m_{\nu}=0.06$ eV
(Table 8). We see that $\sigma(N_{\mathrm{eff}})$ and $\sigma(n_{\mathrm{s}})$
do not change when neutrino mass is fixed in the $\Lambda$CDM +
$N_{\mathrm{eff}}$ model compared to when it is free. In the $\Lambda$CDM
model, $\sigma(n_{\mathrm{s}})$ decreases slightly compared to the models
where $N_{\mathrm{eff}}$ is free. $\sigma(H_{0})$ is about a factor of two
smaller in the $\Lambda$CDM model compared to the $\Lambda$CDM +
$N_{\mathrm{eff}}$ \+ $\sum m_{\nu}$.

In order to see how much improvement is gained from delensing even after
adding DESI BAO, we compare the $1\sigma$ parameter uncertainties obtained
from lensed or delensed CMB-HD $TT$, $TE$, $EE$, and $BB$ power spectra plus
the lensing $\kappa\kappa$ spectrum and DESI BAO data in Table 9. We find that
delensing does improve parameter constraints, with a 20% improvement for
$N_{\mathrm{eff}}$ which is important for pushing the error well below the
0.027 threshold for a spin-0 particle [5].

Finally, in Table 10, we show the forecasted $1\sigma$ parameter constraints
from CMB-HD plus DESI BAO data for different values of the maximum multipole
$\ell_{\mathrm{max}}$ or $L_{\mathrm{max}}$ used in the analysis. At higher
multipoles the data is more sensitive to non-linear effects in the matter
power spectrum, which are more complex to model. Thus we show the parameter
constraints achievable from the linear and semi-linear regimes. We see that to
obtain $\sigma(n_{\mathrm{s}})=0.0013$ and $\sigma(N_{\mathrm{eff}})=0.014$
requires the full range of multipoles out to 20,000.

Parameter | Fiducial | MCMC | Fisher | Fisher / MCMC
---|---|---|---|---
|  |  $1\sigma$ error |  $1\sigma$ error |
$\Omega_{\mathrm{b}}h^{2}$ | $0.022370$ | $0.0000258$ | $0.0000263$ | $1.02$
$\Omega_{\mathrm{c}}h^{2}$ | $0.12000$ | $0.000410$ | $0.000414$ | $1.01$
$\ln(10^{10}A_{\mathrm{s}})$ | $3.044$ | $0.00953$ | $0.00974$ | $1.02$
$n_{\mathrm{s}}$ | $0.9649$ | $0.00133$ | $0.00134$ | $1.01$
$\tau$ | $0.0544$ | $0.00505$ | $0.00515$ | $1.02$
$100\theta_{\mathrm{MC}}$ | $1.040920$ | $0.0000583$ | $0.0000598$ | $1.03$
$N_{\mathrm{eff}}$ | $3.046$ | $0.0139$ | $0.0143$ | $1.03$
$\sum m_{\nu}$ [eV] | $0.060$ | $0.0243$ | $0.0244$ | $1.00$
Table 6: Here we show the $1\sigma$ uncertainties for a $\Lambda$CDM + $N_{\mathrm{eff}}$ \+ $\sum m_{\nu}$ model from delensed CMB-HD $TT$, $TE$, $EE$, $BB$ and $\kappa\kappa$ power spectra plus DESI BAO. The first two columns list the parameter and its fiducial value; the third column lists the marginalized $1\sigma$ uncertainty on each parameter obtained from an MCMC run, and the fourth column lists the uncertainty estimated from a Fisher analysis. In the last column we take the ratio of the two errors, finding good agreement. |  | $\ell_{\mathrm{max}},\leavevmode\nobreak\ L_{\mathrm{max}}=20,000$ | $\ell_{\mathrm{max}},L_{\mathrm{max}}=10,000$
---|---|---|---
Parameter | Fiducial | HD Lensed | HD Lensed | HD Delensed |  HD Delensed | HD Delensed
$\Lambda$CDM |  |  | \+ FG | \+ FG |  \+ FG + DESI BAO | \+ FG + DESI BAO
$\Omega_{\mathrm{b}}h^{2}$ |  $0.022370$ |  $0.000017$ |  $0.000018$ |  $0.000017$ |  $0.000017$ | $0.000017$
$\Omega_{\mathrm{c}}h^{2}$ |  $0.12000$ |  $0.00046$ |  $0.00048$ |  $0.00045$ |  $0.00037$ | $0.00038$
$\ln(10^{10}A_{\mathrm{s}})$ |  $3.0440$ |  $0.0078$ |  $0.0080$ |  $0.0074$ |  $0.0062$ | $0.0063$
$n_{\mathrm{s}}$ |  $0.9649$ |  $0.0013$ |  $0.0014$ |  $0.0013$ |  $0.0012$ | $0.0013$
$\tau$ |  $0.0544$ |  $0.0045$ |  $0.0046$ |  $0.0043$ |  $0.0036$ | $0.0036$
$H_{0}$ [km s-1 Mpc-1] |  $67.36$ |  $0.17$ |  $0.18$ |  $0.17$ |  $0.14$ | $0.14$
Table 7: Forecasted cosmological parameter constraints from CMB-HD $TT$, $TE$, $EE$, $BB$ and $\kappa\kappa$ power spectra for a $\Lambda$CDM model. All forecasts include a $\tau$ prior of $\tau=0.054\pm 0.007$ from Planck [1]. The first two columns list the parameters and their fiducial values. The following two columns list their forecasted marginalized $1\sigma$ uncertainties when using lensed spectra with or without foregrounds in the temperature maps. The fifth column shows the forecast when delensing the CMB spectra, and the sixth column shows the change when including DESI BAO data [16]. In each of these cases we use CMB and CMB lensing multipoles out to $\ell_{\mathrm{max}},L_{\mathrm{max}}=20,000$ for CMB-HD. The last column lists the same information as the sixth column, but instead using a maximum multipole of $\ell_{\mathrm{max}},L_{\mathrm{max}}=10,000$ for both CMB and CMB lensing power spectra. |  | $\ell_{\mathrm{max}},L_{\mathrm{max}}=20,000$ | $\ell_{\mathrm{max}},L_{\mathrm{max}}=10,000$
---|---|---|---
Parameter | Fiducial | HD Lensed | HD Lensed | HD Delensed |  HD Delensed | HD Delensed
$\Lambda$CDM+$N_{\mathrm{eff}}$ |  |  | \+ FG | \+ FG |  \+ FG + DESI BAO | \+ FG + DESI BAO
$\Omega_{\mathrm{b}}h^{2}$ |  $0.022370$ |  $0.000032$ |  $0.000033$ |  $0.000027$ |  $0.000026$ | $0.000026$
$\Omega_{\mathrm{c}}h^{2}$ |  $0.12000$ |  $0.00051$ |  $0.00052$ |  $0.00047$ |  $0.00041$ | $0.00041$
$\ln(10^{10}A_{\mathrm{s}})$ |  $3.0440$ |  $0.0078$ |  $0.0080$ |  $0.0075$ |  $0.0062$ | $0.0063$
$n_{\mathrm{s}}$ |  $0.9649$ |  $0.0013$ |  $0.0015$ |  $0.0015$ |  $0.0013$ | $0.0014$
$\tau$ |  $0.0544$ |  $0.0045$ |  $0.0046$ |  $0.0043$ |  $0.0036$ | $0.0036$
$H_{0}$ [km s-1 Mpc-1] |  $67.36$ |  $0.21$ |  $0.22$ |  $0.21$ |  $0.18$ | $0.19$
$N_{\mathrm{eff}}$ |  $3.046$ |  $0.016$ |  $0.017$ |  $0.015$ |  $0.014$ | $0.015$
Table 8: Forecasted cosmological parameter constraints from CMB-HD $TT$, $TE$, $EE$, $BB$ and $\kappa\kappa$ power spectra on a $\Lambda$CDM + $N_{\rm{eff}}$ model. All forecasts include a $\tau$ prior of $\tau=0.054\pm 0.007$ from Planck [1]. The first two columns list the parameters and their fiducial values. The following two columns list their forecasted marginalized $1\sigma$ uncertainties when using lensed spectra with or without foregrounds in the temperature maps. The fifth column shows the forecast when delensing the CMB spectra, and the sixth column shows the change when including DESI BAO data [16]. In each of these cases we use CMB and CMB lensing multipoles out to $\ell_{\mathrm{max}},L_{\mathrm{max}}=20,000$ for CMB-HD. The last column lists the same information as the sixth column, but instead using a maximum multipole of $\ell_{\mathrm{max}},L_{\mathrm{max}}=10,000$ for both CMB and CMB lensing power spectra.  Parameter | Fiducial | HD Lensed | HD Delensed
---|---|---|---
$\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$ |  | \+ FG + DESI BAO | \+ FG + DESI BAO
$\Omega_{\mathrm{b}}h^{2}$ |  $0.022370$ |  $0.000032$ | $0.000026$
$\Omega_{\mathrm{c}}h^{2}$ |  $0.12000$ |  $0.00045$ | $0.00041$
$\ln(10^{10}A_{\mathrm{s}})$ |  $3.044$ |  $0.010$ | $0.0098$
$n_{\mathrm{s}}$ |  $0.9649$ |  $0.0014$ | $0.0013$
$\tau$ |  $0.0544$ |  $0.0055$ | $0.0052$
$H_{0}$ [km s-1 Mpc-1] |  $67.36$ |  $0.30$ | $0.29$
$N_{\mathrm{eff}}$ |  $3.046$ |  $0.017$ | $0.014$
$\sum m_{\nu}$ [eV] |  $0.06$ |  $0.025$ | $0.025$
Table 9: Here we compare the forecasted $1\sigma$ parameter constraints from lensed (third column) and delensed (fourth column) CMB-HD $TT$, $TE$, $EE$, $BB$ and $\kappa\kappa$ power spectra when including DESI BAO in order to see how much improvement is gained from delensing alone. The parameter names and fiducial values are listed in the first and second columns, respectively. We find that delensing improves parameter constraints even after including DESI BAO, with the most improvement seen for $N_{\mathrm{eff}}$. Parameter |  |  $\ell_{\mathrm{max}},\leavevmode\nobreak\ L_{\mathrm{max}}$ for HD Delensed + FG + DESI BAO
---|---|---
$\Lambda$CDM+$N_{\mathrm{eff}}$+$\sum m_{\nu}$ | Fiducial | 1000 | 3000 | 5000 | 10,000 | 20,000
$\Omega_{\mathrm{b}}h^{2}$ |  $0.022370$ |  $0.00016$ |  $0.000039$ |  $0.000027$ |  $0.000026$ | $0.000026$
$\Omega_{\mathrm{c}}h^{2}$ |  $0.12000$ |  $0.0024$ |  $0.00054$ |  $0.00042$ |  $0.00041$ | $0.00041$
$\ln(10^{10}A_{\mathrm{s}})$ |  $3.044$ |  $0.014$ |  $0.011$ |  $0.011$ |  $0.010$ | $0.0098$
$n_{\mathrm{s}}$ |  $0.9649$ |  $0.0048$ |  $0.0023$ |  $0.0018$ |  $0.0014$ | $0.0013$
$\tau$ |  $0.0544$ |  $0.0065$ |  $0.0060$ |  $0.0057$ |  $0.0054$ | $0.0052$
$H_{0}$ [km s-1 Mpc-1] |  $67.36$ |  $0.80$ |  $0.32$ |  $0.29$ |  $0.29$ | $0.29$
$N_{\mathrm{eff}}$ |  $3.046$ |  $0.15$ |  $0.030$ |  $0.018$ |  $0.015$ | $0.014$
$\sum m_{\nu}$ [eV] |  $0.06$ |  $0.037$ |  $0.029$ |  $0.028$ |  $0.026$ | $0.025$
Table 10: Shown are the forecasted $1\sigma$ uncertainties from delensed CMB-
HD $TT$, $TE$, $EE$, $BB$ and $\kappa\kappa$ spectra plus DESI BAO, when
varying the maximum multipole $\ell_{\mathrm{max}}$ or $L_{\mathrm{max}}$ for
the CMB spectra. The parameter names and fiducial values are listed in the
first and second columns, respectively, while the remaining columns list the
$1\sigma$ parameter errors for the given maximum multipole. We see that to
obtain $\sigma(n_{\mathrm{s}})=0.0013$ and $\sigma(N_{\mathrm{eff}})=0.014$,
for example, requires the full range of multipoles out to 20,000.

## References

  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Astronomy & Astrophysics 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Aiola  _et al._ [2020] S. Aiola _et al._ (ACT), J. Cosmol. Astropart. Phys. 12, 047 (2020), arXiv:2007.07288 [astro-ph.CO] .
  * Balkenhol  _et al._ [2023] L. Balkenhol _et al._ (SPT-3G), Phys. Rev. D 108, 023510 (2023), arXiv:2212.05642 [astro-ph.CO] .
  * Ade  _et al._ [2019] P. Ade _et al._ (Simons Observatory), J. Cosmol. Astropart. Phys. 02, 056 (2019), arXiv:1808.07445 [astro-ph.CO] .
  * Abazajian  _et al._ [2019] K. Abazajian _et al._ , (2019), arXiv:1907.04473 [astro-ph.IM] .
  * Sehgal  _et al._ [2019] N. Sehgal _et al._ , (2019), arXiv:1906.10134 [astro-ph.CO] .
  * Aiola  _et al._ [2022] S. Aiola _et al._ (CMB-HD), (2022), arXiv:2203.05728 [astro-ph.CO] .
  * Han  _et al._ [2021] D. Han _et al._ (ACT), J. Cosmol. Astropart. Phys. 01, 031 (2021), arXiv:2007.14405 [astro-ph.CO] .
  * Millea  _et al._ [2021] M. Millea _et al._ , Astrophys. J. 922, 259 (2021), arXiv:2012.01709 [astro-ph.CO] .
  * Carron  _et al._ [2017] J. Carron, A. Lewis, and A. Challinor, J. Cosmol. Astropart. Phys. 05, 035 (2017), arXiv:1701.01712 [astro-ph.CO] .
  * Green  _et al._ [2017] D. Green, J. Meyers, and A. van Engelen, J. Cosmol. Astropart. Phys. 12, 005 (2017), arXiv:1609.08143 [astro-ph.CO] .
  * Hotinli  _et al._ [2022] S. C. Hotinli, J. Meyers, C. Trendafilova, D. Green, and A. van Engelen, J. Cosmol. Astropart. Phys. 04, 020 (2022), arXiv:2111.15036 [astro-ph.CO] .
  * Chung  _et al._ [2020] E. Chung, S. Foreman, and A. van Engelen, Phys. Rev. D 101, 063534 (2020), [Erratum: Phys.Rev.D 102, 109903 (2020)], arXiv:1910.09565 [astro-ph.CO] .
  * McCarthy  _et al._ [2022] F. McCarthy, J. C. Hill, and M. S. Madhavacheril, Phys. Rev. D 105, 023517 (2022), arXiv:2103.05582 [astro-ph.CO] .
  * Han and Sehgal [2022] D. Han and N. Sehgal, Phys. Rev. D 105, 083516 (2022), arXiv:2112.02109 [astro-ph.CO] .
  * Aghamousa  _et al._ [2016] A. Aghamousa _et al._ (DESI), (2016), arXiv:1611.00036 [astro-ph.IM] .
  * Torrado and Lewis [2021] J. Torrado and A. Lewis, J. Cosmol. Astropart. Phys. 2021, 057 (2021), arXiv:2005.05290 [astro-ph.IM] .
  * Lewis  _et al._ [2000] A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. 538, 473 (2000), arXiv:astro-ph/9911177 .
  * Howlett  _et al._ [2012] C. Howlett, A. Lewis, A. Hall, and A. Challinor, J. Cosmol. Astropart. Phys. 1204, 027 (2012), arXiv:1201.3654 [astro-ph.CO] .
  * Mandal  _et al._ [2022] S. Mandal, N. Sehgal, and T. Namikawa, Phys. Rev. D 105, 063537 (2022), arXiv:2201.02204 [astro-ph.CO] .
  * Mead  _et al._ [2015] A. Mead, J. Peacock, C. Heymans, S. Joudaki, and A. Heavens, Monthly Notices of the Royal Astronomical Society 454, 1958 (2015), arXiv:1505.07833 [astro-ph.CO] .
  * Mead  _et al._ [2016] A. Mead, C. Heymans, L. Lombriser, J. Peacock, O. Steele, and H. Winther, Monthly Notices of the Royal Astronomical Society 459, 1468 (2016), arXiv:1602.02154 [astro-ph.CO] .
  * Zeldovich and Sunyaev [1969] Y. B. Zeldovich and R. A. Sunyaev, Astrophysics and Space Science 4, 301 (1969).
  * Sunyaev and Zeldovich [1970] R. A. Sunyaev and Y. B. Zeldovich, Astrophysics and Space Science 7, 3 (1970).
  * Sunyaev and Zeldovich [1972] R. A. Sunyaev and Y. B. Zeldovich, Comments on Astrophysics and Space Physics 4, 173 (1972).
  * Foreman  _et al._ [2023] S. Foreman, S. C. Hotinli, M. S. Madhavacheril, A. van Engelen, and C. D. Kreisch, Phys. Rev. D 107, 083502 (2023), arXiv:2209.03973 [astro-ph.CO] .
  * Hu [2001] W. Hu, Astrophysical Journal Letters 557, L79 (2001), arXiv:astro-ph/0105424 .
  * Hu and Okamoto [2002] W. Hu and T. Okamoto, Astrophys. J. 574, 566 (2002), arXiv:astro-ph/0111606 .
  * Okamoto and Hu [2003] T. Okamoto and W. Hu, Phys. Rev. D 67, 083002 (2003), arXiv:astro-ph/0301031 .
  * Hu  _et al._ [2007] W. Hu, S. DeDeo, and C. Vale, New J. Phys. 9, 441 (2007), arXiv:astro-ph/0701276 .
  * Kesden  _et al._ [2003a] M. H. Kesden, A. Cooray, and M. Kamionkowski, Phys. Rev. D 67, 123507 (2003a), arXiv:astro-ph/0302536 .
  * Namikawa  _et al._ [2013] T. Namikawa, D. Hanson, and R. Takahashi, Monthly Notices of the Royal Astronomical Society 431, 609 (2013), arXiv:1209.0091 [astro-ph.CO] .
  * Madhavacheril  _et al._ [2020] M. S. Madhavacheril, K. M. Smith, B. D. Sherwin, and S. Naess,  (2020), 10.1088/1475-7516/2021/05/028, arXiv:2011.02475 [astro-ph.CO] .
  * Nguyen  _et al._ [2019] H. N. Nguyen, N. Sehgal, and M. Madhavacheril, Phys. Rev. D 99, 023502 (2019), arXiv:1710.03747 [astro-ph.CO] .
  * Hanson  _et al._ [2011] D. Hanson, A. Challinor, G. Efstathiou, and P. Bielewicz, Phys. Rev. D 83, 043005 (2011), arXiv:1008.4403 [astro-ph.CO] .
  * Horowitz  _et al._ [2019] B. Horowitz, S. Ferraro, and B. D. Sherwin, Mon. Not. Roy. Astron. Soc. 485, 3919 (2019), arXiv:1710.10236 [astro-ph.CO] .
  * Hadzhiyska  _et al._ [2019] B. Hadzhiyska, B. D. Sherwin, M. Madhavacheril, and S. Ferraro, Phys. Rev. D 100, 023547 (2019), arXiv:1905.04217 [astro-ph.CO] .
  * Chan  _et al._ [2023] V. C. Chan, R. Hložek, J. Meyers, and A. van Engelen, (2023), arXiv:2302.13350 [astro-ph.CO] .
  * Seljak and Hirata [2004] U. Seljak and C. M. Hirata, Phys. Rev. D 69, 043005 (2004), arXiv:astro-ph/0310163 .
  * Smith  _et al._ [2012] K. M. Smith, D. Hanson, M. LoVerde, C. M. Hirata, and O. Zahn, J. Cosmol. Astropart. Phys. 2012, 014 (2012), arXiv:1010.0048 [astro-ph.CO] .
  * Simard  _et al._ [2015] G. Simard, D. Hanson, and G. Holder, Astrophys. J. 807, 166 (2015), arXiv:1410.0691 [astro-ph.CO] .
  * Sehgal  _et al._ [2017] N. Sehgal, M. S. Madhavacheril, B. Sherwin, and A. van Engelen, Phys. Rev. D 95, 103512 (2017), arXiv:1612.03898 [astro-ph.CO] .
  * Benoit-Lévy  _et al._ [2012] A. Benoit-Lévy, K. M. Smith, and W. Hu, Phys. Rev. D 86, 123008 (2012), arXiv:1205.0474 [astro-ph.CO] .
  * Peloton  _et al._ [2017] J. Peloton, M. Schmittfull, A. Lewis, J. Carron, and O. Zahn, Phys. Rev. D 95, 043508 (2017), arXiv:1611.01446 [astro-ph.CO] .
  * Schmittfull  _et al._ [2013] M. M. Schmittfull, A. Challinor, D. Hanson, and A. Lewis, Phys. Rev. D 88, 063012 (2013), arXiv:1308.0286 [astro-ph.CO] .
  * Manzotti  _et al._ [2014] A. Manzotti, W. Hu, and A. Benoit-Lévy, Phys. Rev. D 90, 023003 (2014), arXiv:1401.7992 [astro-ph.CO] .
  * Kesden  _et al._ [2003b] M. H. Kesden, A. Cooray, and M. Kamionkowski, Phys. Rev. D 67, 123507 (2003b), arXiv:astro-ph/0302536 .
  * Eisenstein  _et al._ [2005] D. J. Eisenstein _et al._ (SDSS), Astrophys. J. 633, 560 (2005), arXiv:astro-ph/0501171 .
  * Allison  _et al._ [2015] R. Allison, P. Caucal, E. Calabrese, J. Dunkley, and T. Louis, Phys. Rev. D 92, 123535 (2015), arXiv:1509.07471 [astro-ph.CO] .
  * Dodelson [2003] S. Dodelson,  _Modern Cosmology_ (Academic Press, Amsterdam, 2003).
  * Albrecht  _et al._ [2006] A. Albrecht _et al._ , (2006), arXiv:astro-ph/0609591 .
  * Huterer and Takada [2005] D. Huterer and M. Takada, Astropart. Phys. 23, 369 (2005), arXiv:astro-ph/0412142 .
  * Amara and Refregier [2008] A. Amara and A. Refregier, Monthly Notices of the Royal Astronomical Society 391, 228 (2008), arXiv:0710.5171 [astro-ph] .
  * Bernal  _et al._ [2020] J. L. Bernal, N. Bellomo, A. Raccanelli, and L. Verde, J. Cosmol. Astropart. Phys. 10, 017 (2020), arXiv:2005.09666 [astro-ph.CO] .
  * Hogg and Foreman-Mackey [2018] D. W. Hogg and D. Foreman-Mackey, Astrophysical Journal Supplement 236, 11 (2018), arXiv:1710.06068 [astro-ph.IM] .
  * Lewis and Bridle [2002] A. Lewis and S. Bridle, Phys. Rev. D 66, 103511 (2002), arXiv:astro-ph/0205436 .
  * Lewis [2013] A. Lewis, Phys. Rev. D 87, 103529 (2013), arXiv:1304.4473 [astro-ph.CO] .
  * Gelman and Rubin [1992] A. Gelman and D. B. Rubin, Statistical Science 7, 457 (1992).
  * Tröster  _et al._ [2022] T. Tröster _et al._ , Astronomy & Astrophysics 660, A27 (2022), arXiv:2109.04458 [astro-ph.CO] .
  * Mead  _et al._ [2021] A. Mead, S. Brieden, T. Tröster, and C. Heymans, Monthly Notices of the Royal Astronomical Society 502, 1401 (2021), arXiv:2009.01858 [astro-ph.CO] .
  * McCarthy  _et al._ [2021] F. McCarthy, S. Foreman, and A. van Engelen, Phys. Rev. D 103, 103538 (2021), arXiv:2011.06582 [astro-ph.CO] .
  * Wallisch [2018] B. Wallisch,  _Cosmological Probes of Light Relics_ , Ph.D. thesis, Cambridge U. (2018), arXiv:1810.02800 [astro-ph.CO] .
  * Green  _et al._ [2019] D. Green _et al._ , Bull. Am. Astron. Soc. 51, 159 (2019), arXiv:1903.04763 [astro-ph.CO] .
  * Peccei [2008] R. D. Peccei, Lect. Notes Phys. 741, 3 (2008), arXiv:hep-ph/0607268 .
  * Marsh [2018] D. J. E. Marsh, in  _13th Patras Workshop on Axions, WIMPs and WISPs_ (2018) pp. 59–74, arXiv:1712.03018 [hep-ph] .
  * Di Luzio  _et al._ [2020] L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, Phys. Rept. 870, 1 (2020), arXiv:2003.01100 [hep-ph] .
  * Baumann  _et al._ [2016] D. Baumann, D. Green, and B. Wallisch, Phys. Rev. Lett. 117, 171301 (2016), arXiv:1604.08614 [astro-ph.CO] .
  * Blum  _et al._ [2014] K. Blum, R. T. D’Agnolo, M. Lisanti, and B. R. Safdi, Phys. Lett. B 737, 30 (2014), arXiv:1401.6460 [hep-ph] .
  * Hannestad  _et al._ [2007] S. Hannestad, A. Mirizzi, G. G. Raffelt, and Y. Y. Y. Wong, J. Cosmol. Astropart. Phys. 08, 015 (2007), arXiv:0706.4198 [astro-ph] .
  * Di Valentino  _et al._ [2016] E. Di Valentino, E. Giusarma, M. Lattanzi, O. Mena, A. Melchiorri, and J. Silk, Phys. Lett. B 752, 182 (2016), arXiv:1507.08665 [astro-ph.CO] .
  * Schaan and Ferraro [2019] E. Schaan and S. Ferraro, Phys. Rev. Lett. 122, 181301 (2019), arXiv:1804.06403 [astro-ph.CO] .
  * Legrand and Carron [2022] L. Legrand and J. Carron, Phys. Rev. D 105, 123519 (2022), arXiv:2112.05764 [astro-ph.CO] .
  * Sailer  _et al._ [2023] N. Sailer, S. Ferraro, and E. Schaan, Phys. Rev. D 107, 023504 (2023), arXiv:2211.03786 [astro-ph.CO] .
  * Legrand and Carron [2023] L. Legrand and J. Carron, (2023), arXiv:2304.02584 [astro-ph.CO] .