11institutetext: International School for Advanced Studies (SISSA), Via
Bonomea 265, I-34136 Trieste, Italy 22institutetext: National Institute for
Nuclear Physics (INFN) – Sezione di Trieste, Via Valerio 2, I-34127 Trieste,
Italy 33institutetext: Department of Physics, University of Oxford, Denys
Wilkinson Building, Keble Road, Oxford OX1 3RH, United Kingdom
44institutetext: Kavli Institute for the Physics and Mathematics of the
Universe (Kavli IPMU, WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba
277-8583, Japan 55institutetext: Department of Physics, Florida State
University, Tallahassee, Florida 32306, USA 66institutetext: Instituto de
Astrofísica and Centro de Astro-Ingeniería, Facultad de Física, Pontificia
Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul,
Santiago, Chile 77institutetext: Université Paris Cité, CNRS, Astroparticule
et Cosmologie, F-75013 Paris, France 88institutetext: Institute for
Fundamental Physics of the Universe (IFPU), Via Beirut 2, I-34151 Grignano
(TS), Italy 99institutetext: Berkeley Center for Cosmological Physics,
Department of Physics, University of California, Berkeley, CA 94720, USA
1010institutetext: Lawrence Berkeley National Laboratory, One Cyclotron Road,
Berkeley, CA 94720, USA 1111institutetext: Jodrell Bank Centre for
Astrophysics, School of Physics and Astronomy, The University of Manchester,
Oxford Road, Manchester M20 4PE, United Kingdom 1212institutetext: School of
Physics and Astronomy, Cardiff University, The Parade, Cardiff, Wales CF24
3AA, UK 1313institutetext: Center for Computational Astrophysics, Flatiron
Institute, 162 5th Avenue, 10010, New York, NY, USA 1414institutetext:
Department of Physics, University of Texas at Austin, Austin, Texas 78722, USA
1515institutetext: Kavli Institute for the Physics and Mathematics of the
Universe (Kavli IPMU,WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba
277-8583,Japan 1616institutetext: Department of Physics, University of Milano-
Bicocca, Piazza della Scienza 3, 20126 Milan, Italy 1717institutetext: Joseph
Henry Laboratories of Physics, Jadwin Hall, Princeton University, Princeton,
NJ, USA 08544 1818institutetext: Department of Astrophysical Sciences, Peyton
Hall, Princeton University, Princeton, NJ USA 08544

# The Simons Observatory: pipeline comparison and validation for large-scale
$B$-modes

Kevin Wolz  email: K. Wolz<EMAIL_ADDRESS>Susanna Azzoni 3344 Carlos
Hervías-Caimapo 5566 Josquin Errard 77 Nicoletta Krachmalnicoff 112288 David
Alonso 33 Carlo Baccigalupi 112288 Antón Baleato Lizancos 991010 Michael L.
Brown 1111 Erminia Calabrese 1212 Jens Chluba 1111 Jo Dunkley 17171818 Giulio
Fabbian 12121313 Nicholas Galitzki 1414 Baptiste Jost 771515 Magdy Morshed 77
Federico Nati 1616

###### Abstract

Context. The upcoming Simons Observatory Small Aperture Telescopes aim at
achieving a constraint on the primordial tensor-to-scalar ratio $r$ at the
level of $\sigma(r=0)\lesssim 0.003$, observing the polarized CMB in the
presence of partial sky coverage, cosmic variance, inhomogeneous non-white
noise, and Galactic foregrounds.

Aims. We present three different analysis pipelines able to constrain $r$
given the latest available instrument performance, and compare their
predictions on a set of sky simulations that allow us to explore a number of
Galactic foreground models and elements of instrumental noise, relevant for
the Simons Observatory.

Methods. The three pipelines use different combinations of parametric and non-
parametric component separation at the map and power spectrum levels, and
employ $B$-mode purification to estimate the CMB $B$-mode power spectrum. They
are tested and compared regarding their capability to analyze a common set of
simulated realistic frequency maps, and to extract constraints on the tensor-
to-scalar ratio $r$. Their performance is evaluated in terms of bias and
statistical uncertainty on this parameter.

Results. In most of the scenarios the three methodologies achieve similar
performance. Nevertheless, several simulations with complex foreground signals
lead to a $>2\sigma$ bias on $r$ if analyzed with the default versions of
these pipelines, highlighting the need for more sophisticated pipeline
components that marginalize over foreground residuals. We show two such
extensions, using power-spectrum-based and map-based methods, that are able to
fully reduce the bias on $r$ below the statistical uncertainties in all
foreground models explored, at a moderate cost in terms of $\sigma(r)$.

##  1 Introduction

One of the next frontiers in cosmological science using the Cosmic Microwave
Background (CMB) is the observation of large-scale $B$-mode polarization, and
the consequential potential detection of primordial gravitational waves. Such
a detection would let us glance into the early Universe and its very high
energy physics, at scales unattainable by any other experiment. Primordial
tensor perturbations, which would constitute a stochastic background of
primordial gravitational waves, would source a parity-odd $B$-mode component
in the polarization of the CMB (Kamionkowski et al., 1997; Seljak, 1997;
Seljak & Zaldarriaga, 1997; Zaldarriaga & Seljak, 1997). The ratio between the
amplitudes of the primordial power spectrum of these tensor perturbations and
the primordial spectrum of the scalar perturbations is referred to as the
tensor-to-scalar ratio $r$. This ratio covers a broad class of models of the
early Universe, allowing us to test and discriminate between models that
predict a wide range of values of $r$. These include vanishingly small values,
as resulting from models of quantum gravity (e.g. Ijjas & Steinhardt, 2018,
2019), as well as those expected to soon enter the detectable range, predicted
by models of inflation (Starobinskiǐ, 1979; Abbott & Wise, 1984; Martin et
al., 2014b, a; Planck Collaboration X, 2020). An unequivocal measurement of
$r$, or a stringent upper bound, would thus greatly constrain the landscape of
theories on the early Universe.

Although there is no evidence of primordial $B$-modes yet, current CMB
experiments have placed stringent constraints on their amplitude, finding
$r<0.036$ at 95% confidence (BICEP/Keck Collaboration, 2021) when evaluated at
a pivot scale of 0.05 Mpc-1. At the same time, these experiments have firmly
established that the power spectrum of primordial scalar perturbations is not
exactly scale-independent, with the scalar spectral index $n_{s}-1\sim 0.03$
(e.g. Planck Collaboration VI, 2020). Given this measurement, several classes
of inflationary models predict $r$ to be in the $\sim 10^{-3}$ range (see
Kamionkowski & Kovetz, 2016, and references therein).

Even though the only source of primordial large-scale $B$-modes at linear
order are tensor fluctuations, in practice a measurement is complicated by
several factors: first, the gravitational deflection of the background CMB
photons by the cosmic large-scale structure creates coherent sub-degree
distortions in the CMB, known as CMB lensing (Lewis & Challinor, 2006).
Through this mechanism, the nonlinear scalar perturbations from the late
Universe transform a fraction of the parity-even $E$-modes into $B$-modes at
intermediate and small scales (Zaldarriaga & Seljak, 1998). Second, diffuse
Galactic foregrounds have significant polarized emission, and in particular
foreground components such as synchrotron radiation and thermal emission from
dust produce $B$-modes with a significant amplitude. Component separation
methods, which exploit the different spectral energy distributions (SED) of
the CMB and foregrounds to separate the different components, are thus of
vital importance (Delabrouille & Cardoso, 2007; Leach et al., 2008). Practical
implementations of these methods must also be able to carry out this
separation in the presence of instrumental noise and systematic effects (e.g.
Natoli et al., 2018; Abitbol et al., 2021).

Polarized Galactic foregrounds pose a formidable obstacle when attempting to
measure primordial $B$-modes at the level of $r\sim 10^{-3}$. Current
measurements of Galactic emission demonstrate that the Galactic $B$-mode
signal is dominant over the cosmological signal on the relevant scales (Planck
Collaboration X, 2016; Planck Collaboration Int. XXX, 2016; Planck
Collaboration IV, 2020; Planck Collaboration XI, 2020). At the minimum of
polarized Galactic thermal dust and synchrotron, around 80 GHz, their $B$-mode
signal represents an effective tensor-to-scalar ratio with amplitude larger
than the target CMB signal, even in the cleanest regions of the sky
(Krachmalnicoff et al., 2016). Component separation methods are able to clean
most of this, but small residuals left after the cleaning could be comparable
to the primordial $B$-mode signal we want to measure. In recent years, many
works have analyzed this problem and made forecasts of how well we could
potentially measure $r$ with different ground-based and satellite experiments
(e.g. Betoule et al., 2009; Bonaldi & Ricciardi, 2011; Katayama & Komatsu,
2011; Armitage-Caplan et al., 2012; Errard & Stompor, 2012; Remazeilles et
al., 2016; Stompor et al., 2016; Errard et al., 2016; Hervías-Caimapo et al.,
2017; Alonso et al., 2017; Remazeilles et al., 2018b, a; Errard & Stompor,
2019; Thorne et al., 2019; Remazeilles et al., 2021; Azzoni et al., 2021;
Hervías-Caimapo et al., 2022; CMB-S4 Collaboration, 2022; Vacher et al., 2022;
LiteBIRD Collaboration, 2022). In general, these works have highlighted how,
if left untreated, systematic residuals arising from a simplistic
characterization of foregrounds will bias an $r\sim 10^{-3}$ measurement by
several $\sigma$. Thus, it is of vital importance to model the required
foreground complexity when cleaning the multi-frequency CMB observations and
to keep a tight control of systematics without introduction of bias.

Multiple upcoming CMB experiments rank the detection of large-scale primordial
$B$-modes among their primary science targets. Near-future experiments such as
the BICEP Array (Hui et al., 2018) target a detection at the level of $r\sim
0.01$, while in the following decade, next-generation projects such as
LiteBIRD (Hazumi et al., 2019) and CMB-S4 (Abazajian et al., 2016) aim at
$r\sim 0.001$.

The Simons Observatory (SO), like the BICEP Array, targets the detection of
primordial gravitational waves at the level of $r\sim 0.01$ (see “The Simons
Observatory: science goals and forecasts”, SO Collaboration, 2019), and its
performance at realizing this goal is the main focus of this paper. SO is a
ground-based experiment, located at the Cerro Toco site in the Chilean Atacama
desert, which will observe the microwave sky in six frequency channels, from
27 to 280 GHz, due to start in 2023. SO consists of two main instruments. On
the one hand, a Large Aperture Telescope (LAT) with a 6m diameter aperture
will target small-scale CMB physics, secondary anisotropies and the CMB
lensing signal. Measurements of the latter will serve to subtract lensing-
induced $B$-modes from the CMB signal to retrieve primordial $B$-modes (using
a technique known as “delensing”, see Namikawa et al., 2022). On the other
hand, multiple Small Aperture Telescopes (SATs) with 0.5m diameter apertures
will make large-scale, deep observations of $\sim 10$% of the sky, with the
main aim of constraining the primordial $B$-mode signal, peaking on scales
$\ell\sim 80$ (the so-called “recombination bump”). See SO Collaboration
(2019) for an extended discussion on experimental capabilities.

In this paper, we aim at validating three independent $B$-mode analysis
pipelines. We compare their performance regarding a potential $r$ measurement
by the SO SATs, and evaluate the capability of the survey to constrain
$\sigma(r=0)\leq 0.003$ in the presence of foreground contamination and
instrumental noise. To that end, we produce sky simulations encompassing
different levels of foreground complexity, CMB with different values of $r$
and different amounts of residual lensing contamination, and various levels of
the latest available instrumental noise111We note that the pipelines are still
agnostic to some aspects of the instrumental noise such as filtering, which
may impact the overall forecasted scientific performance. We anticipate
studying these in detail in future work., calculated from the parametric
models presented in SO Collaboration (2019).

We feed these simulations through the analysis pipelines and test their
performance, quantifying the bias and statistical uncertainty on $r$ as a
function of foreground and noise complexity. The three pipelines are described
in detail in Section 2. Section 3 presents the simulations used in the
analysis, including the models used to produce CMB and foreground sky maps, as
well as instrumental noise. In Section 4, we present our forecasts for $r$,
accompanied by the power spectrum products and a channel weights comparison.
Section 5 shows preliminary results on a set of new, complex foreground
simulations. In Section 6 we summarize and draw our conclusions. Finally,
Appendix A summarizes the $\chi^{2}$ analysis performed on the
Cross-$C_{\ell}$ cleaning pipeline, while Appendix B discusses biases on
Gaussian simulations observed with the NILC cleaning pipeline.

##  2 Methods, pipelines

In this section we present our three component separation pipelines, that
adopt complementary approaches widely used in the literature: power-spectrum-
based parametric cleaning (BICEP2 Collaboration & Keck Array Collaboration
(2016, 2018)), Needlet Internal Linear Combination (NILC) blind cleaning
(Delabrouille et al., 2009; Basak & Delabrouille, 2012, 2013), and map-based
parametric cleaning (Poletti & Errard, 2023). In the following, these are
denominated pipelines A, B and C, respectively. The cleaning algorithms
operate on different data spaces (harmonic, needlet and pixel space) and vary
in their cleaning strategy (parametric, meaning that we assume an explicit
model for the frequency spectrum of the foreground components, or blind,
meaning that we do not model the foregrounds or make any assumptions on how
their frequency spectrum should be). Hence, they do not share the same set of
method-induced systematic errors. This will serve as an important argument in
favor of claiming robustness of our inference results.

Table 1 lists the three pipelines and their main properties. Although there
are some similarities between these analysis pipelines and the forecasting
frameworks that were exploited in SO Collaboration (2019), the tools developed
for this paper are novel implementations designed to deal with realistic SO
data-like inputs, including complex noise and more exotic foreground
simulations compared to what was considered in the previous work. We stress
again that no filtering or other systematic effects were included in the noise
maps.

Table 1: Overview of the component separation pipelines used to infer $r$ Pipeline | method | data space | blind / parametric |  $r$ inference step
---|---|---|---|---
A | Cross-$C_{\ell}$ cleaning | harmonic (power spectra) | parametric | multi-frequency $C_{\ell}$ likelihood
B | NILC cleaning | needlets (maps) | blind | CMB-only $C_{\ell}$ likelihood
C | map-based cleaning | pixels (maps) | parametric | CMB-only $C_{\ell}$ likelihood

###  2.1 Pipeline A: Cross-$C_{\ell}$ cleaning

Pipeline A is based on a multi-frequency power-spectrum-based component
separation method, similar to that used in the latest analysis carried out by
the BICEP/Keck collaboration (BICEP2 Collaboration & Keck Array Collaboration,
2016, 2018). The data vector is the full set of cross-power spectra between
all frequency maps, $C_{\ell}^{\nu\nu^{\prime}}$. The likelihood compares this
against a theoretical prediction that propagates the map-level sky and
instrument model onto the corresponding power spectra. The full pipeline is
publicly available222See github.com/simonsobs/BBPower., and a schematic
overview is provided in Figure 1.

In step 1, power spectra are measured using a pseudo-$C_{\ell}$ approach with
$B$-mode purification as implemented in NaMaster (Alonso et al., 2019). As
described in Smith & Zaldarriaga (2007) the presence of a sky mask leads to
the presence of ambiguous modes contaminated by full-sky $E$-modes. These must
be removed at the map level to avoid the contribution to the power spectrum
uncertainties from the leaked $E$-modes. The mask used for this analysis
traces the hits count map released in SO Collaboration (2019) (see Figure 2),
and its edges are apodized using a C1-type kernel (see Grain et al., 2009)
with an apodization scale of 10 degrees, yielding an effective sky coverage of
$f_{\rm sky}\sim 10\%$. Each power spectrum is calculated in bandpower windows
with constant bin width $\Delta\ell=10$, of which we only keep the range
$30\leq\ell\leq 300$. Our assumption is that on real data, larger scales are
contaminated by atmospheric noise and filtering, whereas smaller scales,
targeted by the SO-LAT and useful for constraining lensing $B$-modes, do not
contain any significant primordial $B$-mode contribution. To avoid a
significant bias in the auto-correlations when removing the impact of
instrumental noise, a precise noise model is required that may not be
available in practice. We address this issue by using data splits, which in
the case of real data may be formed by subdividing data among observation
periods, sets of detectors, sky patches or by other means, while in this
paper, we resort to simulations.

We construct simulated observations for each sky realization comprising $S=4$
independent splits with the same sky but different noise realizations (each
with a commensurately larger noise amplitude). We compute $BB$ power spectra
from pairs of maps, each associated with a given data split and a given
frequency channel. For any fixed channel pair combination, we average over the
corresponding set of $S(S-1)/2=6$ power spectra with unequal split pairings.
For $N=6$ SAT frequency channels, this results in a collection of
$N(N+1)/2=21$ noise-debiased multi-frequency power spectra, shown in Figure 3.
Note that we could have modeled and subtracted the noise bias explicitly,
since we have full control over the noise properties in our simulations. In
realistic settings, however, the accuracy of an assumed noise model may be
limited. While inaccurate noise modeling would affect the statistical
uncertainty $\sigma(r)$ through the covariance matrix calculated from
simulations, the cross-split approach ensures robustness of the inferred value
of $r$ against noise-induced bias.

In step 2, we estimate the bandpower covariance matrix from simulations,
assuming no correlations between different multipole windows. Note that, to
include foreground signal variance in the budget, Gaussian foreground
simulations (see Section 3) are considered for the covariance computation
step, since our realistic foreground templates cannot be used as statistical
samples. As we show in Appendix A, this covariance matrix is indeed
appropriate, as it leads to the theoretically expected empirical distribution
of the $\hat{\chi}^{2}_{\rm min}$ statistic not only in the case of Gaussian
foregrounds, but also for the non-Gaussian foreground simulations. Inaccurate
covariance estimates would make this statistic peak at higher or lower values,
which we do not observe.

Step $3$ is the parameter inference stage. We use a Gaussian likelihood when
comparing the multi-frequency power spectra with their theoretical prediction.
Note that, in general, the power spectrum likelihood is non-Gaussian, and
Hamimeche & Lewis (2008) provide an approximate likelihood that is able to
account for this non-Gaussianity. We explicitly verified that both likelihoods
lead to equivalent parameter constraints, and thus choose the simpler Gaussian
option. The validity of the Gaussian approximation is a consequence of the
central limit theorem, since each measured bandpower consists of effectively
averaging over $N_{\rm mode}\simeq\Delta\ell\times
f_{\rm{sky}}\times(2\ell+1)>61$ independent squared modes on the scales used
here. Note that this assumption is valid thanks to the relatively large SO-SAT
sky patch and would not hold any longer for BICEP/Keck-like patch sizes. The
default sky model is the same as that described in Abitbol et al. (2021). We
model the angular power spectra of dust and synchrotron as power laws of the
form $D_{\ell}=A_{c}(\ell/\ell_{0})^{\alpha_{c}}$, with $\ell_{0}=80$, and
$c=d$ or $s$ for dust and synchrotron respectively. The dust SED is modeled as
a modified black-body spectrum with spectral index $\beta_{d}$ and temperature
$\Theta_{d}$, which we fix to $\Theta_{d}=20\,{\rm K}$. The synchrotron SED is
modeled as a power law with spectral index $\beta_{s}$. Finally, we consider a
dust-synchrotron correlation parameter $\epsilon_{ds}$. Including the tensor-
to-scalar ratio $r$ and a free lensing $B$-mode amplitude $A_{\rm lens}$, this
fiducial model has $9$ free parameters:

$\\{A_{\rm
lens},r,A_{d},\alpha_{d},\beta_{d},A_{s},\alpha_{s},\beta_{s},\epsilon_{ds}\\}.$
(1)

We will refer to results using this model as “$C_{\ell}$-fiducial”. Table 2
lists the priors on its $9$ parameters.

The main drawback of power-spectrum-based pipelines in their simplest
incarnation, is their inability to account for spatial variation in the
foreground spectra. If ignored, this spatial variation can lead to biases at
the level of $r\sim O(10^{-3})$, which are significant for the SO target. At
the power spectrum level, spatially-varying SEDs give rise to frequency
decorrelation, which can be included in the model. Here, we will show results
for an extended model that uses the moment-based 333See Tegmark (1998); Chluba
et al. (2017); Vacher et al. (2023) for more details on the moment-expansion
formalism in the context of CMB foregrounds, and Mangilli et al. (2021) as an
alternative power-spectrum-based description. parameterization of Azzoni et
al. (2021) to describe the spatial variation of $\beta_{d}$ and $\beta_{s}$.
The model introduces 4 new parameters

$\\{B_{s},\gamma_{s},B_{d},\gamma_{d}\\}\,,$ (2)

where $B_{c}$ parameterizes the amplitude of the spatial variations in the
spectral index of component $c$, and $\gamma_{c}$ is their power spectrum
slope (see Azzoni et al. (2021) for further details). We will refer to results
using this method as “$C_{\ell}$-moments”, or “A + moments”. The priors in the
shared parameter space are the same as in $C_{\ell}$-fiducial, and Table 2
lists the priors on its additional $4$ parameters. For both methods, we sample
posteriors using the emcee code (Foreman-Mackey et al., 2013). Note that we
assume a top-hat prior on $r$ in the range $[-0.1,\,0.1]$ instead of imposing
$r>0$. The reason is that we would like to remain sensitive to potential
negative biases on $r$. While negative $r$ values do not make sense
physically, they may result from e.g. volume effects caused by choosing
specific priors on other parameters that we have marginalized over. Opening
the prior on $r$ to negative values allows us to monitor these unwanted
effects, offering a simple robustness check. On real data, this will be
replaced by a positivity prior $r>0$, but only after ensuring that our
specific prior choices on the other parameters do not bias $r$, which is the
focus of a future work.

Table 2: Parameter priors for pipeline A, considering both the $C_{\ell}$-fiducial model and the $C_{\ell}$-moments model. Prior types are either Gaussian (G) or top-hat (TH), considered distributed symmetrically around the center value with half width meaning the standard deviation (Gaussian) or the half width (top-hat).  model |  $C_{\ell}$-fiducial and $C_{\ell}$-moments |  $C_{\ell}$-moments only
---|---|---
parameter | $A_{\rm lens}$ | $r$ | $A_{d}$ | $\alpha_{d}$ | $\beta_{d}$ | $A_{s}$ | $\alpha_{s}$ | $\beta_{s}$ | $\epsilon_{ds}$ | $B_{s}$ | $\gamma_{s}$ | $B_{d}$ | $\gamma_{d}$
prior type | TH | TH | TH | TH | G | TH | TH | G | TH | TH | TH | TH | TH
center value | $1.0$ | $0.0$ | $25$ | $0.0$ | $1.54$ | $2.0$ | $-1.0$ | $-3.0$ | $0.0$ | $0.0$ | $-4.0$ | $5.0$ | $-4.0$
half width | $1.0$ | $0.1$ | $25$ | $0.5$ | $0.11$ | $2.0$ | $1.0$ | $0.3$ | $1.0$ | $10.0$ | $2.0$ | $5.0$ | $2.0$
Figure 1: Schematic of pipeline A. Orange colors mark steps that are repeated
500 times, once for each simulation. Figure 2: Apodized SAT hits map with
effective $f_{\rm sky}=10\%$ used in this paper, shown in equatorial
projection. Its edges are apodized using a C1-type kernel with an apodization
scale of 10 degrees.  Figure 3: Power spectrum data analyzed by pipeline A,
showing a single realization of CMB signal and Gaussian foregrounds. Blue
shaded areas quantify the 1$\sigma$ Gaussian standard deviation calculated
from simulations of CMB, noise and Gaussian foregrounds. Note that negative
auto-spectra can occur at noise-dominated scales as a result of cross-
correlating data splits.

###  2.2 Pipeline B: NILC cleaning

Our second pipeline is based on the blind Internal Linear Combination (ILC)
method, which assumes no information on foregrounds whatsoever, and instead
only assumes that the observed data contains one signal of interest (the CMB)
plus noise and contaminants (Bennett et al., 2003). The method assumes a
simple model for the observed multi-frequency maps $\mathbf{d}_{\nu}$ at
$N_{\nu}$ frequency channels (in either pixel or harmonic space)

$\mathbf{d}_{\nu}=a_{\nu}\mathbf{s}+\mathbf{n}_{\nu}\text{,}$ (3)

where $a_{\nu}$ is the black-body spectrum of the CMB, $\mathbf{s}$ is the
amplitude of the true CMB signal and $\mathbf{n}_{\nu}$ is the contamination
in channel $\nu$, which includes the foregrounds and instrumental noise. ILC
exploits the difference between the black-body spectrum of the CMB and the
SED(s) of other components that may be present in the data. The method aims at
reconstructing a map of the CMB component $\tilde{\mathbf{s}}$ as a linear
combination of the data with a set of weights $w_{\nu}$ allowed to vary across
the map,

$\tilde{\mathbf{s}}=\sum_{\nu}w_{\nu}\mathbf{d}_{\nu}=\mathbf{w}^{T}\hat{\sf
d}\,,$ (4)

where both $\mathbf{w}$ and $\hat{\sf d}$ are $N_{\nu}\times N_{\rm pix}$
matrices,with $N_{\rm pix}$ being the number of pixels. Optimal weights are
found by minimizing the variance of $\tilde{\bf s}$. The result is given by:

$\mathbf{w}^{T}=\frac{\mathbf{a}^{T}\hat{\sf C}^{-1}}{\mathbf{a}^{T}\hat{\sf
C}^{-1}\mathbf{a}}\text{,}$ (5)

where $\mathbf{a}$ is the black-body spectrum of the CMB (i.e. a vector filled
with ones if maps are in thermodynamic temperature units) and $\hat{\sf
C}=\langle\hat{\sf d}\,\hat{\sf d}^{T}\rangle$ is the frequency-frequency
covariance matrix per pixel of the observed data. Note that we assume no
correlation between pixels for the optimal weights.

In our particular implementation, we use the Needlet Internal Linear
Combination method (NILC, Delabrouille et al., 2009; Basak & Delabrouille,
2012, 2013). NILC uses localization in pixel and harmonic space by finding
different weights ${\bf w}$ for a set of harmonic filters, called “needlet
windows”. These windows are defined in harmonic space $h_{i}(\ell)$ for
$i=0,...,n_{\rm windows}-1$ and must satisfy the constraint
$\sum_{i=0}^{n_{\rm windows}-1}h_{i}(\ell)^{2}=1$ in order to preserve the
power of the reconstructed CMB. We use $n_{\rm windows}=5$ needlet windows
shown in Figure 4, and defined by

$h_{i}(\ell)=\begin{cases}\cos(\frac{\pi}{2}(\ell^{\rm
peak}_{i}-\ell)/(\ell^{\rm peak}_{i}-\ell^{\rm min}_{i}))&\text{if $\ell^{\rm
min}_{i}\leq\ell<\ell^{\rm peak}_{i}$ }\\\ 1&\text{if $\ell=\ell^{\rm
peak}_{i}$}\\\ \cos(\frac{\pi}{2}(\ell-\ell^{\rm peak}_{i})/(\ell^{\rm
max}_{i}-\ell^{\rm peak}_{i}))&\text{if $\ell^{\rm peak}_{i}<\ell\leq\ell^{\rm
max}_{i}$}\end{cases}\text{,}$ (6)

with $\ell_{\rm min}=\\{$0, 0, 100, 200, 350$\\}$, $\ell_{\rm max}=\\{$100,
200, 350, 500, 500$\\}$, and $\ell_{\rm peak}=\\{$0, 100, 200, 350, 500$\\}$
for the corresponding 5 needlet windows. Even though we do not use the full
500-$\ell$ range where windows are defined for the likelihood sampling, we
still perform the component separation on all 5 windows up to multipoles
beyond our upper limit of $\ell=300$, in order to avoid edge effects on the
smaller scales.

Let us now describe the NILC procedure as illustrated in Figure 5. In step 1,
we perform our CMB reconstruction in the $E$ and $B$ field instead of $Q$ and
$U$. We transform the observed maps to $a_{\ell m}^{X}$ with $X\in E,B$. All
frequency channels are then brought to a common beam resolution by rescaling
the harmonic coefficients with an appropriate harmonic beam window function.
The common beam we adopt is the one from the third frequency channel at 93
GHz, which corresponds to a FWHM of 30 arcmin. For each needlet window index
$i$, we multiply $a_{\ell m}^{X}$ with $h_{i}(\ell)$ as a harmonic filter.
Since different frequency channels have different limiting resolutions, we do
not use all channels in every needlet window: the first 2 windows use all 6
frequency channels, the third window does not use the 27 GHz channel and the
last 2 needlet windows do not use the 27 and 39 GHz channels. The covariance
matrix $\hat{\sf C}$ has dimensions $N_{\nu}\times N_{\nu}\times N_{\rm pix}$.
For each pixel $p$, its corresponding $N_{\nu}\times N_{\nu}$ elements are
computed directly from the data, averaging over the pixels inside a given
pixel domain $\mathcal{D}(p,i)$ around each pixel. In practice, the element
$\nu,\nu^{\prime}$ of the covariance matrix is calculated by multiplying the
two filtered maps at channels $\nu$ and $\nu^{\prime}$, then smoothing that
map with a Gaussian kernel with FWHM equal to the size of the pixel domain
$\mathcal{D}(p,i)$. The FWHMs for the pixel domain size are 185, 72, 44, 31,
and 39 degrees for each needlet window, respectively. 444The domain sizes are
estimated directly from the needlet window scale (see details in the appendix
A of Delabrouille et al., 2009). The ILC bias can be minimized by enlarging
the pixel domains to be big enough to include a higher number of modes. We
choose the resulting ILC bias to not exceed $0.2\%$, for which we need pixel
domain sizes large enough so that each needlet window contains at least 2500
modes.

We then proceed to calculate the weights $\mathbf{w}^{T}$ (see Equation 5) for
window $i$, which is an array with shape (2, $N_{\nu}$, $N_{\rm pixels}^{i}$),
with the first dimension corresponding to the $E$ and $B$ fields. Note that
the number of pixels $N_{\rm pixels}^{i}$ is different for each needlet
window, since we use different pixel resolutions depending on the smallest
scale covered by the window. Finally, we apply Equation 4 to obtain an ILC-
reconstructed CMB map for window $i$. The final step is to filter this map in
harmonic space for a second time with the $h_{i}(\ell)$ window. The final
reconstructed CMB map is the sum of these maps for all five needlet windows.

Figure 4: The five needlet windows used in pipeline B. The black dashed lines
are the beam transfer functions $b_{\ell}^{\nu}$ for the six SO-SAT frequency
channels. The FWHM of the beams is listed in Table 3. The gray shaded area are
the multipole ranges we do not use in the likelihood sampling to estimate $r$.
Figure 5: Schematic of pipeline B. Orange colors mark steps that are repeated
500 times, once for each simulation.

In step 2, the reconstructed CMB maps are compressed into power spectra using
NaMaster, deconvolving to the common beam resolution. We use $B$-mode
purification as implemented in the software and the mask shown in Figure 2. We
estimate the noise bias $N_{\ell}$ in the final map by computing the power
spectrum of noise-only simulations processed with the needlet weights and
windows obtained from the simulated data as described above. $N_{\ell}$ is
averaged over simulations and subtracted from the $C_{\ell}$ of the
reconstructed maps.

Finally, in step 3, we run a Monte Carlo Markov Chain (MCMC) over the
reconstructed $BB$ spectrum (we ignore $EE$ and $EB$) with only two free
parameters: the tensor-to-scalar ratio $r$ and the amplitude of the $BB$
lensing spectrum $A_{\rm lens}$, using the Python package emcee. Both
parameters have a top hat prior (between 0 and 2 for $A_{\rm lens}$, and
between $-0.013$ and infinity for $r$). The covariance matrix is calculated
directly over 500 simulations with the same setup but with Gaussian
foregrounds. As likelihood, we use the same Gaussian likelihood used in
pipeline A and restrict the inference to a multipole range $30<\ell\leq 300$.

While the NILC implementation described above is blind, it can be extended to
a semi-blind approach that introduces a certain level of foreground modeling.
For example, constrained ILC (cILC, Remazeilles et al., 2011) explicitly
nullifies one or more contaminant components in the observed data (such as
thermal dust) by including their modeled SED in the variance minimization
problem that calculates the weights in Equation 5. This foreground modeling
can be further extended to include the moment expansion of the SED described
in Section 2.1. This method, known as constrained moment ILC (cMILC
Remazeilles et al., 2021), has been shown to be effective for cleaning the
large-scale $B$-mode contamination for space experiments such as LiteBIRD.
While not used in this work, these extensions and others will be considered in
future analyses with more complex foregrounds and systematics.

###  2.3 Pipeline C: map-based cleaning

Figure 6: Schematic of pipeline C. Orange colors indicate repetition for each
simulation.

Our third pipeline is a map-based parametric pipeline based on the fgbuster
code Poletti & Errard (2023). This approach is based on the following data
model:

$\displaystyle\centering{\bf d}=\hat{\sf A}{\bf s}+{\bf n}\@add@centering$ (7)

where ${\bf d}$ is a vector containing the polarized frequency maps, $\hat{\sf
A}=\hat{\sf A}(\boldsymbol{\beta})$ is the so-called mixing matrix assumed to
be parameterized by a set of spectral indices $\boldsymbol{\beta}$, ${\bf s}$
is a vector containing the $Q$ and $U$ amplitudes of the sky signals (CMB,
foregrounds) and ${\bf n}$ is the noise contained in each frequency map.
Starting from the input observed data sets ${\bf d}$, Figure 6 shows a
schematic of the pipeline, which contains four steps.

Step 0 is the preprocessing of input simulations: for each simulation, we
combine the simulated noise maps, the foreground and CMB maps and save them on
disk. We create a new set of frequency maps, $\mathbf{\bar{d}}$, smoothed with
a common Gaussian kernel of $100\arcmin$ FWHM.

Step 1 is the actual component separation stage. We optimize the spectral
likelihood defined as Stompor et al. (2009):

$-2\log\left(\mathcal{L}_{\rm spec}(\boldsymbol{\beta})\right)=\left(\hat{\sf
A}^{T}\hat{\sf N}^{-1}\bar{\bf d}\right)^{T}\left(\hat{\sf A}^{T}\hat{\sf
N}^{-1}\hat{\sf A}\right)^{-1}\left(\hat{\sf A}^{T}\hat{\sf N}^{-1}\bar{\bf
d}\right)$ (8)

which uses the common resolution frequency maps, $\bar{\bf d}$, built during
step 0. The right hand side of Equation 8 contains a sum over the observed sky
pixels, assumed to have uncorrelated noise – the diagonal noise covariance
matrix $\hat{\sf N}$ is computed from 500 noise-only simulations. Although, in
principle, $\hat{\sf N}$ can be non-diagonal, we do not observe any
significant bias of the spectral likelihood due to this approximation in this
study. By minimizing Equation 8 we estimate the best-fit spectral indices
$\tilde{\boldsymbol{\beta}}$ and the corresponding mixing matrix $\tilde{\sf
A}\equiv\hat{\sf A}(\tilde{\boldsymbol{\beta}})$. We also estimate the
uncertainties on the recovered spectral indices as provided by the minimizer,
a truncated-Newton algorithm (Nash, 1984) as implemented in scipy (Virtanen et
al., 2020). Having thus obtained estimators of the foreground SEDs, we can
recover the sky component maps with the generalized least-square equation

$\tilde{\bf s}=\left(\tilde{\sf A}^{T}\hat{\sf N}^{-1}\tilde{\sf
A}\right)^{-1}\tilde{\sf A}^{T}\hat{\sf N}^{-1}{\bf d}\equiv\hat{\sf W}{\bf
d}\;,$ (9)

where ${\bf d}$ is the input raw data, and not the common resolution maps. In
steps 1 and 2, we have the possibility to use an inhomogeneous noise
covariance matrix i.e. $\hat{\sf N}=\hat{\sf N}(\hat{\bf n})$ and, although
this is not exploited in this work, a spatially varying mixing matrix
$\boldsymbol{\beta}=\boldsymbol{\beta}(\hat{\bf n})$. For the latter, one can
use the multi-patch or clustering methods implemented in fgbuster (Errard &
Stompor, 2019; Puglisi et al., 2022).

Step 2 comprises the calculation of angular power spectra. The recovered CMB
polarization map is transformed to harmonic space using NaMaster. We estimate
an effective transfer function, ${\bf B}^{\rm eff}_{\ell}=\hat{\sf W}{\bf
B}_{\ell}$ associated with the reconstructed components $\tilde{\bf s}$, from
the channel-specific beams ${\bf B}_{\ell}$. Correcting for the impact of this
effective beam is vital to obtain an unbiased $BB$ spectrum of the foreground-
cleaned CMB, $\tilde{C}_{\ell}^{\rm CMB}$. In the second step, we also
estimate the noise bias from noise simulations, i.e.

$\tilde{\sf N}_{\ell}=\frac{1}{\rm N_{sim}}\sum_{\rm
sims}\sum_{m=-\ell}^{\ell}\frac{\tilde{\bf n}_{\ell,m}\tilde{\bf
n}_{\ell,m^{\prime}}^{\dagger}}{2\ell+1}.$ (10)

where $\tilde{\bf n}=\hat{\sf W}{\bf n}^{\rm sim}$ is the noise in the
recovered component-separated sky maps. We consider 500 simulations to
estimate the noise bias.

Step 3 is the cosmological analysis stage. We model the angular power spectrum
of the component-separated CMB map, including the noise contribution, as

$C_{\ell}^{\rm CMB}(r,A_{\rm lens})\equiv C_{\ell}^{\rm prim}(r)+C_{\ell}^{\rm
lens}(A_{\rm lens})+\tilde{N}^{\rm CMB}_{\ell}$ (11)

and compare data and model with the cosmological likelihood

$\displaystyle-2\log{\mathcal{L}^{\rm cosmo}}$
$\displaystyle\qquad=\sum_{\ell}\left(2\ell+1\right)f_{\rm
sky}\left(\frac{\tilde{C}_{\ell}^{\rm CMB}}{C_{\ell}^{\rm
CMB}}+\log(C_{\ell}^{\rm CMB})\right).$ (12)

It is worth noting that this is only an approximation to the true map-level
Gaussian likelihood which approximates the effective number of modes in each
multipole after masking and purification as $f_{\rm sky}(2\ell+1)$, thus
neglecting any mode-coupling effects induced by the survey footprint. We grid
the likelihood above along the two dimensions $r$ and $A_{\rm lens}$. For each
simulation we then estimate the maximum-likelihood values and 68% credible
intervals from the marginal distributions of $r$ and $A_{\rm lens}$. We
verified that the distributions of recovered $\\{r,\,A_{\rm lens}\\}$ across
simulations are well described by a Gaussian, hence supporting the Gaussian
likelihood in Equation 12.

Pipeline C also offers the option to marginalize over a dust template. The
recovered components in $\tilde{\bf s}$, Equation 9, include the dust $Q$ and
$U$ maps which are typically recovered with high signal-to-noise. In the same
way that we compute $\tilde{C}_{\ell}^{\rm CMB}$ in step 2, we compute the
$BB$ component of the recovered dust map, $\tilde{C}_{\ell}^{\rm dust}$. We
then update our cosmological likelihood, Equation 11, by adding a dust term:

$C_{\ell}^{\rm CMB}=C_{\ell}^{\rm CMB}(r,\,A_{\rm lens})+A_{\rm
dust}\tilde{C}_{\ell}^{\rm dust}.$ (13)

This is a similar approach to earlier methods (Errard & Stompor, 2019;
LiteBIRD Collaboration, 2022). When choosing this approach, the inference of
$r$ during step 3 therefore involves the marginalization over both parameters
$A_{\rm lens}$ and $A_{\rm dust}$. In principle one could add synchrotron or
other terms in Equation 13 but we limit ourselves to dust as it turns out to
be the largest contamination, and, in practice, marginalizing over it allows
us to get unbiased estimates of cosmological parameters. In the remainder of
this paper, we will will refer to this method as “C + dust marginalization”.

##  3 Description of input simulations

Table 3: Instrument and noise specifications used to produce the simulations in this work. Note that these levels correspond to homogeneous noise, while our default analysis assumes noise maps weighted according to the SAT hits map. |  | Baseline | Goal | Pessimistic | Optimistic |
---|---|---|---|---|---|---
Frequency [GHz] | FWHM [arcmin] |  $N_{\rm white}$ [$\mu$K-arcmin] |  $N_{\rm white}$ [$\mu$K-arcmin] | $\ell_{\rm knee}$ | $\ell_{\rm knee}$ | $\alpha_{\rm knee}$
27 | 91 | 46 | 33 | 30 | 15 | -2.4
39 | 63 | 28 | 22 | 30 | 15 | -2.4
93 | 30 | 3.5 | 2.5 | 50 | 25 | -2.5
145 | 17 | 4.4 | 2.8 | 50 | 25 | -3.0
225 | 11 | 8.4 | 5.5 | 70 | 35 | -3.0
280 | 9 | 21 | 14 | 100 | 40 | -3.0

We build a set of dedicated simulations on which to test our data analysis
pipelines and compare results. The simulated maps include cosmological CMB
signal, Galactic foreground emission as well as instrumental noise.

###  3.1 Instrumental specifications and noise

We simulate polarized Stokes $Q$ and $U$ sky maps as observed by the SO-SAT
telescopes. All maps are simulated using the HEALPix pixelation scheme (Górski
et al., 2005) with resolution parameter $N_{\rm side}=512$.

We model the SO-SAT noise power spectra as

$N_{\ell}=N_{\rm white}\Big{[}1+\Big{(}\frac{\ell}{\ell_{\rm
knee}}\Big{)}^{\alpha_{\rm knee}}\Big{]},$ (14)

where $N_{\rm white}$ is the white noise component while $\ell_{\rm knee}$,
and $\alpha_{\rm knee}$ describe the contribution from $1/f$ noise. Following
SO Collaboration (2019) (hereinafter SO2019), we consider different possible
scenarios: “baseline” and “goal” levels for the white noise component, and
“pessimistic” and “optimistic” correlated noise. The values of $(N_{\rm
white},\ell_{\rm knee},\alpha_{\rm knee})$ associated with the different cases
are reported in Table 3. Note that the values of $N_{\rm white}$ correspond to
noise on a sky fraction of $f_{\rm sky}=10\%$ and 5 years of observation time,
as in SO2019. Differently from SO2019, we cite polarization noise levels
corresponding to a uniform map coverage, accounting for the factor of $\sim
1.3$ difference compared to Table 1 in SO2019. We simulate noise maps as
Gaussian realizations of the $N_{\ell}$ power spectra. In our main analysis,
we use noise maps with pixel weights computed from the SO-SAT hits map (see
Figure 2) and refer to this as “inhomogeneous noise”. In Section 4.1, we
briefly present results obtained from equally weighted noise pixels, which we
refer to as “homogeneous noise”. Otherwise, all results in this paper assume
inhomogeneous noise. Note that, although inhomogeneous, the noise realizations
used here lack some of the important anisotropic properties of realistic $1/f$
noise, such as stripes due to the scanning strategy. Thus, together with the
impact of other time-domain effects (e.g. filtering), we leave a more thorough
study of the impact of instrumental noise properties for future work.

###  3.2 CMB

We simulate the CMB signal as Gaussian random realizations following a power
spectrum given by the Planck 2018 best-fit $\Lambda$CDM parameters. Our
baseline model does not include any primordial tensor signal ($r=0$) but
incorporates lensing power in the $BB$ spectra ($A_{\rm lens}=1$). We consider
also two modifications of this model: (i) primordial tensor signal with
$r=0.01$, representing a $\gtrsim 3\sigma$ target detection for SO with
$\sigma(r)=0.003$, as forecasted by SO2019; (ii) reduced lensing power with
$A_{\rm lens}=0.5$, corresponding to a 50% delensing efficiency, achievable
for SO (Namikawa et al., 2022).

For all the different scenarios we simulate 500 realizations of the CMB
signal, convolved with Gaussian beams for each frequency channel, with FWHMs
as reported in Table 3.

###  3.3 Foregrounds

Thermal emission from Galactic dust grains and synchrotron radiation are known
to be the two main contaminants to CMB observations in polarization, at
intermediate and large angular scales, impacting therefore measurements of the
primordial $BB$ signal. Many studies have been performed in the past years on
the characterization of polarized Galactic foreground emission, thanks to the
analysis of WMAP and Planck data, as well as low frequency surveys (Harper et
al., 2022; Krachmalnicoff et al., 2018). However, many aspects of their
emission remain unconstrained including, in particular, the characterization
of their SEDs and its variation across the sky. For this reason, to properly
assess the impact of foreground emission on component separation and $r$
constraints we use four sets of models to simulate the sky emission. As
specified in the following, we use the Python sky model (PySM) package (Thorne
et al., 2017) to simulate polarized foreground components, with some
additional modifications:

  * •

Gaussian foregrounds: we simulate thermal dust emission and synchrotron
radiation as Gaussian realizations of power law $EE$ and $BB$ power spectra.
Although inaccurate, since foregrounds are highly non-Gaussian, this idealized
model has been used to validate the different pipelines and to build
approximate signal covariance matrices from 500 random realizations. In
particular, we estimate the amplitudes of the foreground signal (evaluated for
$D_{\ell}=\ell(\ell+1)C_{\ell}/2\pi$ at $\ell=80$) and the slope of their
polarization angular power spectra from the synchrotron and thermal dust
templates of PySM, evaluated in the SO-SAT sky patch, leading to the following
values ($d$: thermal dust at 353 GHz; $s$: synchrotron at 23 GHz):
$A^{d}_{EE}=56$ $\mu K_{\rm CMB}^{2}$, $A^{d}_{BB}=28$ $\mu K_{\rm CMB}^{2}$,
$\alpha^{d}_{EE}=-0.32$, $\alpha^{d}_{BB}=-0.16$; $A^{s}_{EE}=9$ $\mu K_{\rm
CMB}^{2}$, $A^{s}_{BB}=1.6$ $\mu K_{\rm CMB}^{2}$, $\alpha^{s}_{EE}=-0.7$,
$\alpha^{s}_{BB}=-0.93$. The frequency scaling of the maps at the SO
frequencies has been considered to be a modified black body for thermal dust
emission, with fixed spectral parameter $\beta_{d}=1.54$ and $T_{d}=20$ K, and
a power law for synchrotron with $\beta_{s}=-3$ (in antenna temperature
units).

  * •

d0s0 model: in this case multi-frequency maps have been obtained using the
d0s0 PySM model. This model includes templates for thermal dust emission
coming from Planck high frequency observations and from WMAP 23 GHz maps from
synchrotron radiation. SEDs are considered to be uniform across the sky with
the same values of the spectral parameters used for the Gaussian simulations.

  * •

d1s1 model: this model uses the same foreground amplitude templates as d0s0,
but with the inclusion of spatial variability for spectral parameters, as
described in Thorne et al. (2017).

  * •

dmsm model: this model represents a modification of the d1s1 spatial variation
of spectral parameters. For thermal dust we smoothed the $\beta_{d}$ and
$T_{d}$ templates at an angular resolution of 2 degrees, in order to down-
weight the contribution of instrumental noise fluctuations in the original
PySM maps. For synchrotron emission we modified the $\beta_{s}$ PySM in order
to account for the additional information coming from the analysis of S-PASS
data at 2.3 GHz (see Krachmalnicoff et al. (2018)). In particular S-PASS data
show that the synchrotron spectral index presents enhanced variations with
respect to the PySM template. We therefore multiplied the fluctuations in the
$\beta_{s}$ map by a factor 1.6 to take into consideration larger variations.
Moreover, small scale fluctuations (with a minimum angular resolution of 2
degrees) have been added, as Gaussian realization of a power-law power
spectrum with slope $-2.6$ (see Figure 11 in Krachmalnicoff et al. (2018)).

Note that this quartet of foreground models includes the model d1s1, used for
large-scale $B$-mode forecasts in SO Collaboration (2019), and therefore
represents an extension of the Galactic foreground scenarios probed in that
past analysis. As for the CMB simulations, the multi-frequency foreground maps
at the SO reference frequencies have been convolved with Gaussian beams, to
reach the expected angular resolution. We assume delta-like frequency
bandpasses in order to accelerate the production of these simulations,
although all pipelines are able to handle finite bandpasses. Therefore this
approximation should not impact the performance of any of the pipelines
presented here.

##  4 Results and discussion

Simulations were generated for 4 different noise and foreground models,
respectively, (see Section 3), for a total of 16 different foreground-noise
combinations. For the main analysis, we consider a fiducial CMB model with
$A_{\rm lens}=1$ (no delensing) and $r=0$ (no primordial tensor fluctuations).
In addition, we explored three departures from the fiducial CMB model, with
input parameters $(A_{\rm lens}=0.5,\,r=0)$, $(A_{\rm lens}=1,\,r=0.01)$ and
$(A_{\rm lens}=0.5,\,r=0.01)$. Here we report the results found for all these
cases.

###  4.1 Constraints on $r$

Figure 7: Compilation of the mean $r$ with its $(16,\,84)$% credible interval
as derived from 500 simulations, applying the three nominal pipelines (plus
extensions) to four foreground scenarios of increasing complexity. We assume a
fiducial cosmology with $r=0$ and $A_{\rm lens}=1$, inhomogeneous noise with
goal sensitivity and optimistic $1/f$ noise component (dot markers), and
inhomogeneous noise with baseline sensitivity and pessimistic $1/f$ noise
component (cross markers). Note that the NILC results for Gaussian foregrounds
are based on a smaller sky mask, see Appendix B. Table 4: Results on the mean
of $r$ and the $(16,\,84)$% credible interval derived from 500 simulations, as
inferred by three pipelines (and two extensions) on four foreground models and
four noise cases, two of which are shown in Figure 7. No delensing is assumed
in these fiducial results.

$10^{3}\times(r\pm\sigma(r))$
---
Noise | FG model | Pipeline A |  $+$ moments | Pipeline B | Pipeline C |  $+$ dust marg.
| Gaussian | $-0.1\pm 2.1$ | $0.0\pm 2.8$ | $0.6\pm 2.6^{\dagger}$ | $1.6\pm 2.7$ | $-1.8\pm 4.4$
Goal rms, | d0s0 | $-0.4\pm 2.1$ | $-0.5\pm 2.7$ | $-0.1\pm 2.1$ | $0.5\pm 2.3$ | $-1.7\pm 3.5$
optimistic $1/f$ | d1s1 | $1.8\pm 2.1$ | $-0.2\pm 2.8$ | $2.1\pm 2.1$ | $2.6\pm 2.3$ | $0.0\pm 3.3$
| dmsm | $3.9\pm 2.1$ | $0.3\pm 2.7$ | $3.8\pm 2.1$ | $5.3\pm 2.4$ | $0.2\pm 3.0$
| Gaussian | $-0.2\pm 2.5$ | $-0.1\pm 2.7$ | $1.1\pm 3.3^{\dagger}$ | $0.9\pm 2.1$ | $-0.9\pm 5.3$
Goal rms, | d0s0 | $-0.6\pm 2.5$ | $-0.5\pm 2.8$ | $-0.5\pm 2.8$ | $0.1\pm 2.5$ | $-0.9\pm 4.0$
pessimistic $1/f$ | d1s1 | $1.3\pm 2.5$ | $0.1\pm 3.0$ | $1.2\pm 2.8$ | $3.4\pm 3.1$ | $-0.0\pm 3.9$
| dmsm | $3.2\pm 2.6$ | $0.3\pm 3.9$ | $2.1\pm 2.8$ | $5.5\pm 2.4$ | $0.6\pm 4.2$
| Gaussian | $-0.1\pm 2.6$ | $-0.3\pm 3.3$ | $0.5\pm 3.3^{\dagger}$ | $0.5\pm 3.2$ | $-1.9\pm 5.9$
Baseline rms, | d0s0 | $-0.4\pm 2.6$ | $-0.3\pm 3.3$ | $-0.9\pm 2.7$ | $0.7\pm 2.9$ | $-1.8\pm 4.4$
optimistic $1/f$ | d1s1 | $1.7\pm 2.6$ | $-0.2\pm 3.4$ | $1.0\pm 2.7$ | $1.8\pm 2.7$ | $-0.8\pm 4.8$
| dmsm | $3.9\pm 2.6$ | $0.3\pm 3.5$ | $2.5\pm 2.7$ | $5.5\pm 3.2$ | $0.4\pm 5.0$
| Gaussian | $-0.3\pm 3.4$ | $-0.3\pm 3.8$ | $1.6\pm 4.1^{\dagger}$ | $0.0\pm 2.9$ | $-1.6\pm 5.3$
Baseline rms, | d0s0 | $-0.7\pm 3.4$ | $-0.06\pm 3.9$ | $-0.6\pm 3.6$ | $1.1\pm 3.2$ | $-1.1\pm 5.3$
pessimistic $1/f$ | d1s1 | $1.1\pm 3.4$ | $-0.6\pm 4.0$ | $0.5\pm 3.6$ | $3.8\pm 3.2$ | $-1.2\pm 5.3$
| dmsm | $2.8\pm 3.4$ | $-0.6\pm 4.0$ | $1.2\pm 3.6$ | $6.0\pm 3.1$ | $-0.5\pm 5.1$

† These results are calculated on a smaller, more homogeneous mask, shown in
Figure 14. This is explained in Appendix B.

Let us start by examining the final constraints on $r$ obtained by each
pipeline applied to 500 simulations. These results are summarized in Figure 7
and Table 4. Figure 7 shows the mean $r$ and $(16,\,84)$% credible intervals
found by each pipeline as a function of the input foreground model (labels on
the $x$ axis). Results are shown for 5 pipeline setups: pipeline A using the
$C_{\ell}$-fiducial model (red), pipeline A using the $C_{\ell}$-moments model
(yellow), pipeline B (blue), pipeline C (green) and pipeline C including the
marginalization over the dust amplitude parameter (cyan). For each pipeline,
we show two points with error bars. The dot markers and smaller error bars
correspond to the results found in the best-case instrument scenario (goal
noise level, optimistic $1/f$ component), while the cross markers and larger
error bars correspond to the baseline noise level and pessimistic $1/f$
component. The quantitative results are reported in Table 4.

We will first discuss the nominal pipelines A, B and C without considering any
extensions. We find that for the simpler Gaussian and d0s0 foregrounds, the
nominal pipelines obtain unbiased results, as expected. Pipeline B shows a
slight positive bias for Gaussian foregrounds, in combination with
inhomogeneous noise only. This bias is absent for homogeneous noise and can be
traced back to the pixel covariance matrix used to construct the NILC weights.
We will discuss this in more detail in Appendix B. For now, we show the
results using a more constrained mask that is closer to homogeneity. We stress
that these results, marked with a †, are not comparable to the rest in Table
4, since they are calculated on a different mask. The more complex d1s1
foregrounds lead to a $\sim 1\sigma$ bias in the goal and optimistic noise
scenario. The dmsm foregrounds lead to a noticeable increase of the bias of up
to $\sim 2\sigma$, seen with pipeline C in all noise scenarios and with
pipeline A in the goal-optimistic case, and slightly less with pipeline B. The
modifications introduced in the dmsm foreground model include a larger spatial
variation in the synchrotron spectral index $\beta_{s}$ with respect to d1s1,
and are a plausible reason for the increased bias on $r$.

Remarkably, we find that, in their simplest incarnation, all pipelines achieve
comparable statistical uncertainty on $r$, ranging from $\sigma(r)\simeq
2.1\times 10^{-3}$ to $\sigma(r)\simeq 3.6\times 10^{-3}$ (a $70\%$ increase),
depending on the noise model. Changing between the goal and baseline white
noise levels results in an increase of $\sigma(r)$ of $\sim 20-30\%$. Changing
between the optimistic and pessimistic $1/f$ noise has a similar effect on the
results from pipelines A and B, although $\sigma(r)$ does not increase by more
than $10\%$ when changing to pessimistic $1/f$ noise for pipeline C. These
results are in reasonable agreement with the forecasts presented in SO
Collaboration (2019).

Now let us have a look at the pipeline extensions, A + moments and C + dust
marginalization. Notably, in all noise and foreground scenarios, the two
extensions are able to reduce the bias on $r$ to below $1\sigma$. For the
Gaussian and d0s0 foregrounds, we consistently observe a small negative bias
(at the $\sim 0.1\sigma$ level for A + moments and $<0.5\sigma$ for C + dust
marginalization). This bias may be caused by the introduction of extra
parameters that are prior dominated, like the dust template’s amplitude in the
absence of residual dust contamination, or the moment parameters in the
absence of varying spectral indices of foregrounds. If those extra parameters
are weakly degenerate with the tensor-to-scalar ratio, the marginal $r$
posterior will shift according to the choice of the prior on the extra
parameters. The observed shifts in the tensor-to-scalar ratio and their
possible relation with these volume effects will be investigated in a future
work. For the more complex d1s1 and dmsm, both pipeline extensions effectively
remove the bias observed in the nominal pipelines, achieving a $\sim
0.5\sigma$ bias and lower.

The statistical uncertainty $\sigma(r)$ increases for both pipeline
extensions, although by largely different factors. While C + dust
marginalization yields $\sigma(r)$ between $3.0\times 10^{-3}$ and $5.9\times
10^{-3}$, the loss in precision for A + moments is significantly smaller, with
$\sigma(r)$ varying between $2.7\times 10^{-3}$ and $4.0\times 10^{-3}$
depending on the noise scenario, an average increase of $\sim 25\%$ compared
to pipeline A. In any case, within the assumptions made regarding the SO noise
properties, it should be possible to detect a primordial $B$-mode signal with
$r=0.01$ at the 2–3$\sigma$ level with no delensing. The impact of other
effects, such as time domain filtering or anisotropic noise may affect these
forecasts, and will be studied in more detail in the future.

This analysis was repeated for input CMB maps generated assuming $r=0$ or
$0.01$ and $A_{\rm lens}=0.5$ or $1$. For simplicity, in these cases we
considered only the baseline white noise level with optimistic $1/f$ noise,
and the moderately complex d1s1 foreground model. The results can be found in
Figure 8 and Table 5. A 50% delensing efficiency results in a reduction in the
final $\sigma(r)$ by 25-30% for pipelines A and B, $\sim$10-20% for A +
moments, and 0-33% for C + dust marginalization. The presence of primordial
$B$-modes with a detectable amplitude increases the contribution from cosmic
variance to the error budget, with $\sigma(r)$ growing by up to $40\%$ if
$r=0.01$, in agreement with theoretical expectations. Using C + dust
marginalization, and considering no delensing, we even find $\sigma(r)$
decreasing, hinting at the possible breaking of the degeneracy between $r$ and
$A_{\rm dust}$. We conclude that all pipelines are able to detect the $r=0.01$
signal at the level of $\sim 3\sigma$. As before, we observe a 0.5-1.2$\sigma$
bias on the recovered $r$ that is eliminated by both the moment expansion
method and the dust marginalization method.

Figure 8: Results on $r$ and the $(16,\,84)$% credible interval derived from
500 simulations, applying the three nominal pipelines plus extensions,
assuming input models including primordial $B$-modes and 50% delensing
efficiency. We assume baseline noise level with optimistic $1/f$ component and
the d1s1 foreground template. Table 5: Results on $r$ with its $(16,\,84)$%
credible interval derived from 500 simulations, applying the three nominal
pipelines with extensions, assuming input models including primordial
$B$-modes and 50% delensing efficiency. We assume baseline noise level with
optimistic $1/f$ component and d1s1 foregrounds, see Figure 8.
$10^{3}\times(r\pm\sigma(r))$
---
Input CMB | Pipeline A |  $+$ moments | Pipeline B | Pipeline C |  $+$ dust marg.
$(r=0,$ $A_{\rm{lens}}=1)$ | $1.7\pm 2.6$ | $-0.2\pm 3.4$ | $1.0\pm 2.7$ | $1.8\pm 2.7$ | $-0.8\pm 4.8$
$(r=0,$ $A_{\rm{lens}}=0.5)$ | $1.7\pm 2.0$ | $0.7\pm 2.7$ | $0.7\pm 2.4$ | $3.3\pm 2.7$ | $-0.3\pm 3.6$
$(r=0.01,$ $A_{\rm{lens}}=1)$ | $11.5\pm 3.1$ | $10.6\pm 3.2$ | $11.0\pm 3.2$ | $13.0\pm 3.1$ | $13.4\pm 4.3$
$(r=0.01,$ $A_{\rm{lens}}=0.5)$ | $11.8\pm 2.4$ | $10.3\pm 2.8$ | $10.5\pm 3.0$ | $13.0\pm 3.5$ | $9.4\pm 4.0$

Finally, we have explored how cosmological constraints and the pipelines’
performances are affected by noise inhomogeneity resulting from weighting the
noise pixels according to the SO-SAT hits map. The geographical location of
SO, and the size of the SAT field of view, place constraints on the possible
scanning strategies. In particular, it forces SO to target a patch with a
relatively large sky fraction $f_{\rm sky}\sim 0.15$, and surrounded by a
$\sim 10$ degree wide boundary with significantly higher noise (see hits map
in Figure 2). The lower panel of Figure 9 shows the ratio between the values
of $\sigma(r)$ found using inhomogeneous noise realizations and those with
homogeneous noise in the baseline-optimistic noise model with d0s0
foregrounds, averaged over 500 simulations. We see that for all pipeline
scenarios, $\sigma(r)$ increases by $\sim 30\%$ due to the noise
inhomogeneity.

Figure 9: Results on $r$ with $(16,\,84)$% credible levels derived from 500
simulations, applying the three nominal pipelines plus extensions, assuming
the d0s0 foreground scenario with baseline white noise level and optimistic
$1/f$ component. Cross markers with smaller error bars correspond to
homogeneous noise across the SAT field of view, dot markers with larger error
bars correspond to inhomogeneous noise, and the relative increase in
$\sigma(r)$ between both is shown in the bottom panel.

###  4.2 Power spectra

Having presented the constraints on $r$ found by each pipeline, let us now
examine the CMB power spectrum products. Pipelines B and C produce CMB-only
maps and base their inference of $r$ on the resulting power spectra, whereas
pipeline A works directly with the cross-frequency power spectra of the
original multi-frequency maps. Nevertheless, CMB power spectra are an
important data product that every pipeline should be able to provide.
Following the methods presented in Dunkley et al. (2013) and Planck
Collaboration XI (2016); Planck Collaboration V (2020), we use a modified
version of pipeline A that retrieves CMB-only bandpowers from multi-frequency
power spectra, marginalizing over foregrounds with an MCMC sampler as
presented in Section 2.1. Note that this method, originally developed for
high-$\ell$ CMB science, is applicable since we are in the Gaussian likelihood
regime. By re-inserting this cleaned CMB spectrum into a Gaussian likelihood
with parameters $(r,\,A_{\rm lens})$, we obtain constraints that are
consistent with the results shown in Table 4.

Figure 10: CMB-only power spectra resulting from component separation with
pipelines A, B and C showing the non-Gaussian foregrounds scenarios d0s0 (_top
panel_), d1s1 (_middle panel_) and dmsm (_bottom panel_) while considering the
goal-optimistic noise scenario. The different colored markers with error bars
show the mean of 500 simulations and the scatter between them (corresponding
to the statistical uncertainties of a single realization). The dotted lines in
the corresponding colors indicate the best-fit power spectrum model. In the
dmsm case, we show the extended pipeline results from A + moments and C + dust
marginalization with the best-fit models shown as dot-dashed lines. The black
solid line is the input CMB model containing lensing $B$-modes only. Note that
pipeline C only considers multipoles up to $\ell=180$ in the power spectrum
likelihood.

Figure 10 shows the CMB power spectra for the three complex foreground
simulations d0s0, d1s1 and dmsm (upper, middle and lower panel, respectively)
while considering the goal-optimistic noise scenario. The various markers with
error bars denote the measured CMB power spectra and their $1\sigma$ standard
deviation across 500 simulations, while the black solid line denotes the input
CMB power spectrum. Results are shown in gold triangles, blue circles,
turquoise diamonds for pipeline A, B and C respectively. The dotted lines show
the best-fit CMB model for the three nominal pipelines (using the same color
scheme). Only in the dmsm foreground scenario, which is the most complex
considered here, we also show the results from pipeline C + dust
marginalization (dark red squares with error bars) and the best-fit CMB power
spectrum from A + moments (pink dot-dashed line) and C + dust marginalization
(dark red dot-dashed line).

For the nominal pipelines (A, B and C) without extensions, the measured power
spectra display a deviation from the input CMB at low multipoles, increasing
with rising foreground complexity. For dmsm at multipoles $\lesssim 50$, this
bias amounts to about 1.5$\sigma$ and it goes down to less than 0.5$\sigma$ at
$80\lesssim\ell\lesssim 250$. The three pipelines agree reasonably well, while
pipeline A appears slightly less biased for the lowest multipoles. Pipelines B
and C show an additional mild excess of power in their highest multipole bins,
with a $<0.3\sigma$ increase in pipeline C for $130\lesssim\ell\lesssim 170$
and up to 1$\sigma$ for the highest multipole ($\ell=297$) in pipeline B. This
might indicate power leakage from the multiple operations on map resolutions
implemented in pipelines B and C. In pipeline B, these systematics could come
from first deconvolving the multi-frequency maps and then convolving them with
a common beam in order to bring them to a common resolution, whereas in
pipeline C, the leakage is likely due to the linear combination of the multi-
resolution frequency maps following Equation (9). Other multipole powers lie
within the 1$\sigma$ standard deviation from simulations for all three
pipelines.

Both extensions, A + moments and C + dust marginalization, lead to an unbiased
CMB power spectrum model, as shown by the pink and dark red dot-dashed lines
and the square markers in the lower panel of Figure 10. In the case of
pipelines B and C, comparing the best-fit models obtained from the measured
power spectra to the input CMB model, we find sub-sigma bias for all bins with
$\ell>100$. We show, however, that the ability to marginalize over additional
foreground residuals (e.g. the dust-template marginalization in pipeline C) is
able to reduce this bias on all scales, at the cost of increased
uncertainties. Implementing this capability in the blind NILC pipeline B would
likely allow it to reduce the bias observed in the figure.

The SO-SATs are expected to constrain the amplitude of CMB lensing $B$-modes
to an unprecedented precision. As can be seen from Figure 10, individual,
cleaned CMB bandpowers without delensing at multipoles $\ell\gtrsim 150$
achieve a signal-to-noise ratio of about 10, accounting for a combined
precision on the lensing amplitude of $\sigma(A_{\rm lens})\lesssim 0.03$ when
considering multipoles up to $\ell_{\rm max}=300$. This is consistent with the
inference results obtained by pipelines A and B.

###  4.3 Channel weights

Figure 11: Channel weights for pipelines A, B and C. We show the SMICA weights
for 27 different $\ell$-bins calculated from raw, noisy $C_{\ell}$s (pipeline
A, upper panel), pixel-averaged NILC weights for 5 needlet windows (pipeline
B, middle panel) and pixel-averaged weights from parametric map-based
component separation (pipeline C, lower panel). Weights are averaged over
$100$ simulations, shown are goal white $+$ optimistic $1/f$ noise (dashed
lines) as well as baseline white $+$ pessimistic $1/f$ noise (solid lines).
The semitransparent gray areas represent the channel weights’ 1-$\sigma$
standard deviation across $100$ simulations, covering baseline $+$ pessimistic
noise.

Our three baseline pipelines differ fundamentally in how they separate the sky
components. One common feature among all pipelines is the use of six frequency
channels to distinguish components by means of their different SEDs. In Figure
11 we visualize the channel weights as a function of the band center
frequency, showing the pipelines in three vertically stacked panels. In the
upper panel, we show the effective weights for the CMB applied to the noise-
debiased raw power spectra used by pipeline A, distinguishing between weights
for each harmonic bin:

$\mathbf{w_{\ell}}^{T}=\frac{\mathbf{a}^{T}\hat{\sf
C}_{\ell}^{-1}}{\mathbf{a}^{T}\hat{\sf C}_{\ell}^{-1}\mathbf{a}}.$ (15)

Here, $\hat{\sf C}_{\ell}$ is the $6\times 6$ matrix of raw cross-frequency
power spectra from noisy sky maps, $\mathbf{a}$ is a vector of length $6$
filled with ones. This is equivalent to the weights employed by the SMICA
component separation method (Cardoso et al., 2008) and by ILC as explained in
Section 2.2. The middle panel shows the pixel-averaged NILC weights for the
five needlet windows (Figure 4) used in pipeline B. In the lower panel we show
the CMB weights calculated with the map-based component separation, pipeline
C, averaged over the observed pixels to yield an array of $6$ numbers. All
channel weights are averaged over 100 simulations containing CMB ($r=0$,
$A_{\rm lens}=1$), d1s1 foregrounds and two noise models: goal white noise
with optimistic $1/f$ noise is shown as dashed lines, whereas baseline white
noise with pessimistic $1/f$ noise is shown as solid lines. Moreover, the gray
shaded areas quantify the 1-$\sigma$ uncertainty region of these weights in
the baseline/pessimistic case estimated from $100$ simulations. We see from
Figure 11 that the average channel weights agree well between pipelines A, B
and C. Mid-frequency channels at $93$ and $145$ GHz are assigned positive CMB
weights throughout all pipelines, while high- and low-frequency channels tend
to be suppressed owing to a larger dust and synchrotron contamination,
respectively. In more detail, the $280$ GHz channel is given negative weight
in all pipelines, while average weights at $27$, $39$ and $225$ GHz are
negative with pipeline C and either positive or negative in pipelines A and B,
depending on the angular scale. The CMB channel weight tends to consistently
increase for pipelines A and B as a function of multipole, a fact well
exemplified by NILC at $225$ GHz, matching the expectation that the CMB
lensing signal becomes more important at high $\ell$. Overall, Figure 11
illustrates that foregrounds at low and high frequencies are consistently
subtracted by the three component separation pipelines, with the expected
scale dependence in pipelines A and B. Moreover, at every frequency, the
channel weights are non-negligible and of similar size across the pipelines,
meaning that all channels give a relevant contribution to component separation
for all pipelines.

##  5 More complex foregrounds: d10s5

During the completion of this paper, the new PySM3 Galactic foreground models
were made publicly available. In particular, in these new models, templates
for polarized thermal dust and synchrotron radiation have been updated
including the following changes:

  1. 1.

Large-scale thermal dust emission is based on the GNILC maps (Planck
Collaboration IV, 2020), which present a lower contamination from CIB emission
with respect to the d1 model, based on Commander templates.

  2. 2.

For both thermal dust and synchrotron radiation small scales structures are
added by modifying the logarithm of the polarization fraction tensor555See
pysm3.readthedocs.io/en/latest/preprocess-templates..

  3. 3.

Thermal dust spectral parameters are based on GNILC products, with larger
variation of $\beta_{d}$ and $T_{d}$ parameter at low resolution compared to
the d1 model. Small-scale structure is also added as Gaussian realizations of
power-law power spectra.

  4. 4.

The new template for $\beta_{s}$ includes information from the analysis of
S-PASS data (Krachmalnicoff et al., 2018), in a similar way as the one of the
sm model adopted in this work. In addition, small-scale structures are present
at sub-degree angular scales.

These modifications are encoded in the models called d10 and s5 in the updated
version of PySM. Although these models are still to be considered preliminary,
both in terms of their implementation details in PySM666The PySM library is
currently under development, with beta versions including minor modifications
of the foreground templates being realized regularly. In this part of our
analysis we make use of PySM v3.4.0b3.  and in general, being based on
datasets that may not fully involve the unknown level of foreground
complexity, we decided to dedicate an extra section to their analysis. For
computational speed, we ran the five pipeline set-ups on a reduced set of 100
simulations containing the new d10s5 foregrounds template, CMB with a standard
cosmology ($r=0,\,A_{\rm lens}=1$) and inhomogeneous noise in the goal-
optimistic scenario. The resulting marginalized posterior mean and
$(16,\,84)$% credible intervals on $r$, averaged over 100 simulations, are:

$\displaystyle r=0.0194\pm 0.0021$ $\displaystyle(\text{pipeline A})$
$\displaystyle r=0.0025\pm 0.0031$ $\displaystyle(\text{A + moments})$
$\displaystyle r=0.0144\pm 0.0023$ $\displaystyle(\text{pipeline B})$ (16)
$\displaystyle r=0.0220\pm 0.0026$ $\displaystyle(\text{pipeline C})$
$\displaystyle r=-0.0015\pm 0.0051$ $\displaystyle(\text{C + dust marg.})$

Note that the respective bias obtained with pipelines A, B and C are at 9, 6
and 8$\sigma$, quadrupling the bias of the dmsm foreground model. Crucially,
this bias is reduced to $\sim 0.8\sigma$ with the A + moments pipeline, with a
$40\%$ increase in $\sigma(r)$ compared to pipeline A, and $\sim 0.3\sigma$
with the C + dust-marginalization pipeline, with a $95\%$ increase in
$\sigma(r)$ compared to pipeline C. This makes A + moments the unbiased method
with the lowest statistical error.

The $C_{\ell}$-fiducial model achieves minimum $\chi^{2}$ values of $\approx
600\pm 40$. Although this is an increase of $\Delta\chi^{2}\sim 30$ with
respect to the less complex foreground simulations (see Appendix A), the
associated probability to exceed is $p\sim 0.16$ (assuming our null
distribution is a $\chi^{2}$ with $N_{\rm data}-N_{\rm parameters}=567$
degrees of freedom), and therefore it would not be possible to identify the
presence of a foreground bias by virtue of the model providing a bad fit to
the data. The $\chi^{2}$ value we find also confirms that the covariance
matrix calculated from Gaussian simulations is still appropriate for the non-
Gaussian d10s5 template. On the other hand, A + moments achieves best-fit
$\chi^{2}$ values of $\approx 531\pm 33$, with a $<5\%$ decrease in
$\Delta\chi^{2}$ with respect to less complex foreground simulations,
indicating an improved fitting accuracy.

As shown in Figure 12, the relative model odds between $C_{\ell}$-fiducial and
$C_{\ell}$-moments (see Appendix A for more details) vary between $10^{-26}$
and $10^{-2}$, clearly favoring $C_{\ell}$-moments. Out of 100 d10s5
simulations, 99 yield model odds below 1% and 78 below $10^{-5}$. As opposed
to the less complex foreground simulations (d0s0, d1s1 and dmsm), d10s5 gives
strong preference to using the moment expansion in the power spectrum model.
Note that the AIC-based model odds are computed from the differences of
$\chi^{2}$ values that stem from the same simulation seed, while the
$\chi^{2}$ analysis only considers the models individually, which explains why
the former makes for a much more powerful model comparison test.

These results consider only the most optimistic noise scenario. Other cases
would likely lead to larger uncertainty and, as a consequence, lower relative
biases. In this regard, it is highly encouraging to see two pipeline
extensions being able to robustly separate the cosmological signal from
Galactic synchrotron and dust emission with this high-level complexity. This
highlights the importance of accounting for and marginalizing over residual
foreground contamination due to frequency decorrelation for the level of
sensitivity that SO and other next-generation observatories will achieve.

The contrast between the results obtained on the dmsm and d10s5 simulations
gives us an opportunity to reflect on the strategy one should follow when
determining the fiducial component separation method to use in primordial
$B$-mode searches. Although the dmsm model leads to a $\sim 2\sigma$ bias on
$r$ under the simplest component separation algorithms, simple model-selection
metrics are not able to provide significant evidence that a more sophisticated
modeling of foregrounds is needed. The situation changes with d10s5. A
conservative approach is therefore to select the level of complexity needed
for component separation by ensuring that unbiased constraints are obtained
for all existing foreground models consistent with currently available data.
The analysis methods passing this test can then form the basis for the
fiducial $B$-mode constraints. Alternative results can then be obtained with
less conservative component separation techniques, but their goodness of fit
(or any similar model selection metric) should be compared with that of the
fiducial methods. These results should also be accompanied by a comprehensive
set of robustness tests able to identify signatures of foreground
contamination in the data. This will form the basis of a future work. In a
follow-up paper, we will also explore the new set of complex PySM3 foreground
templates in more detail.

##  6 Conclusions

Figure 12: Empirical distribution of the AIC-based relative model odds (see
Appendix A for details) between the $C_{\ell}$-fiducial and the
$C_{\ell}$-moments model, evaluated for 100 simulations and generated with
five different Galactic foreground templates, including the new PySM
foreground model d10s5. We find strong preference for the $C_{\ell}$-moments
model in the d10s5 foreground scenario, and only then.

In this paper, we have presented three different component separation
pipelines designed to place constraints on the amplitude of cosmological
$B$-modes on polarized maps of the SO Small Aperture Telescopes. The pipelines
are based on multi-frequency $C_{\ell}$ parametric cleaning (Pipeline A),
blind Needlet ILC cleaning (Pipeline B), and map-based parametric cleaning
(Pipeline C). We have also studied extensions of pipelines A and C that
marginalize over additional residual foreground contamination, using a moment
expansion or a dust power spectrum template, respectively. We have tested and
compared their performance on a set of simulated maps containing lensing
$B$-modes with different scenarios of instrumental noise and Galactic
foreground complexity. The presence of additional instrumental complexity,
such as time-domain filtering, or anisotropic noise, are likely to affect our
results. The impact of these effects will be more thoroughly studied in future
work.

We find the inferred uncertainty on the tensor-to-scalar ratio $\sigma(r)$ to
be compatible between the three pipelines. While the simpler foreground
scenarios (Gaussian, d0s0) do not bias $r$, spectral index variations can
cause an increased bias of 1–2$\sigma$ if left untreated, as seen with more
complex foreground scenarios (d1s1, dmsm). Modeling and marginalizing over the
spectral residuals is vital to obtain unbiased $B$-mode estimates. The
extensions to pipelines A and C are thus able to yield unbiased estimates on
all foreground scenarios, albeit with a respective increase in $\sigma(r)$ by
$\sim 20\%$ (A + moments) and $>30\%$ (C + dust marginalization). These
results are in good agreement with the forecasts presented in SO Collaboration
(2019).

After testing on simulations with an $r=0.01$ cosmology, we conclude that
under realistic conditions and if the forecasted map-noise levels and
characteristics are achieved, SO should be able to detect a $r=0.01$ signal at
$\sim 2$–$3\sigma$ after 5 years of observation. Inhomogeneous noise from the
SAT map-making scanning strategy brings about $30\%$ increase in $\sigma(r)$
as compared to homogeneous noise. Analyzing the per-channel weights for our
pipelines, we find all frequency channels to be relevant for the CMB signal
extraction and all pipelines to be in good agreement. These forecasts cover
the nominal SO survey, and can be considered pessimistic in the light of
prospective additional SATs that will further improve the sensitivity on large
angular scales.

We have also carried out a preliminary analysis of new, more complex,
foreground models recently implemented in PySM3, in particular the d10s5
foreground template. The much higher level of spatial SED variation allowed by
this model leads to a drastic increase in the bias on $r$ by up to $\sim
9\sigma$, when analyzed with the nominal pipelines A, B and C. Fortunately,
this bias can be reduced to below $1\sigma$ when using A + moments and c +
dust marginalization. These extensions lead to a 40% and 95% degradation of
the error bars, respectively. Our results highlight the importance of
marginalizing over residuals caused by frequency decorrelation for SO-like
sensitivities. Although we have not analyzed d10s5 to the same depth as the
other foreground models presented here, it is encouraging to confirm that we
have the tools at hand to obtain robust, unbiased constraints on the tensor-
to-scalar ratio in the presence of such complex Galactic foregrounds.

In preparation for the data collected by SO in the near future, we will
continue our investigations into Galactic foreground models with other levels
of complexity as the field progresses. Nevertheless, the current work has
shown that the analysis pipelines in place for SO are able to obtain robust
constraints on the amplitude of primordial $B$ modes in the presence of
Galactic foregrounds covering the full range of complexity envisaged by
current, state-of-the-art models.

###### Acknowledgements.

The authors would like to thank Ken Ganga and Arthur Kosowsky for useful
feedback. The group at SISSA acknowledges support from the COSMOS Network of
the Italian Space Agency and the InDark Initiative of the National Institute
for Nuclear Physics (INFN). KW is funded by a SISSA PhD fellowship. SA is
funded by a Kavli/IPMU doctoral studentship. CHC acknowledges NSF award
1815887 and FONDECYT Postdoc fellowship 3220255. DA is supported by the
Science and Technology Facilities Council through an Ernest Rutherford
Fellowship, grant reference ST/P004474. This work is part of a project that
has received funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program (PI: Josquin
Errard, Grant agreement No. 101044073). ABL is a BCCP fellow at UC Berkeley
and Lawrence Berkeley National Laboratory. MLB acknowledges funding from UKRI
and STFC (Grant awards ST/X006344/1 and ST/X006336/1). EC acknowledges support
from the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation programme (Grant agreement No. 849169). JC was
furthermore supported by the ERC Consolidator Grant CMBSPEC (No. 725456) and
the Royal Society as a Royal Society University Research Fellow at the
University of Manchester, UK (No. URF/R/191023). GF acknowledges the support
of the European Research Council under the Marie Skłodowska Curie actions
through the Individual Global Fellowship No. 892401 PiCOGAMBAS. We acknowledge
the use of CAMB (Lewis et al., 2000), healpy (Zonca et al., 2019), numpy
(Harris et al., 2020), matplotlib (Hunter, 2007), emcee (Foreman-Mackey et
al., 2013) and fgbuster (Errard & Stompor, 2019; Puglisi et al., 2022)
software packages.

## References

  * Abazajian et al. (2016) Abazajian, K. N., Adshead, P., Ahmed, Z., et al. 2016, arXiv e-prints, arXiv:1610.02743
  * Abbott & Wise (1984) Abbott, L. F. & Wise, M. B. 1984, Nuclear Physics B, 244, 541
  * Abitbol et al. (2021) Abitbol, M. H., Alonso, D., Simon, S. M., et al. 2021, J. Cosmology Astropart. Phys., 2021, 032
  * Akaike (1974) Akaike, H. 1974, IEEE Transactions on Automatic Control, 19, 716
  * Alonso et al. (2017) Alonso, D., Dunkley, J., Thorne, B., & Næss, S. 2017, Phys. Rev. D, 95, 043504
  * Alonso et al. (2019) Alonso, D., Sanchez, J., Slosar, A., & LSST Dark Energy Science Collaboration. 2019, MNRAS, 484, 4127
  * Armitage-Caplan et al. (2012) Armitage-Caplan, C., Dunkley, J., Eriksen, H. K., & Dickinson, C. 2012, MNRAS, 424, 1914
  * Azzoni et al. (2021) Azzoni, S., Abitbol, M. H., Alonso, D., et al. 2021, J. Cosmology Astropart. Phys., 2021, 047
  * Basak & Delabrouille (2012) Basak, S. & Delabrouille, J. 2012, MNRAS, 419, 1163
  * Basak & Delabrouille (2013) Basak, S. & Delabrouille, J. 2013, MNRAS, 435, 18
  * Bennett et al. (2003) Bennett, C. L., Hill, R. S., Hinshaw, G., et al. 2003, ApJS, 148, 97
  * Betoule et al. (2009) Betoule, M., Pierpaoli, E., Delabrouille, J., Le Jeune, M., & Cardoso, J. F. 2009, A&A, 503, 691
  * BICEP2 Collaboration & Keck Array Collaboration (2016) BICEP2 Collaboration & Keck Array Collaboration. 2016, Phys. Rev. Lett., 116, 031302
  * BICEP2 Collaboration & Keck Array Collaboration (2018) BICEP2 Collaboration & Keck Array Collaboration. 2018, Phys. Rev. Lett., 121, 221301
  * BICEP/Keck Collaboration (2021) BICEP/Keck Collaboration. 2021, Phys. Rev. Lett., 127, 151301
  * Bonaldi & Ricciardi (2011) Bonaldi, A. & Ricciardi, S. 2011, MNRAS, 414, 615
  * Cardoso et al. (2008) Cardoso, J.-F., Martin, M., Delabrouille, J., Betoule, M., & Patanchon, G. 2008, arXiv e-prints, arXiv:0803.1814
  * Chluba et al. (2017) Chluba, J., Hill, J. C., & Abitbol, M. H. 2017, MNRAS, 472, 1195
  * CMB-S4 Collaboration (2022) CMB-S4 Collaboration. 2022, ApJ, 926, 54
  * Delabrouille & Cardoso (2007) Delabrouille, J. & Cardoso, J. F. 2007, in International Summer School on Data Analysis in Cosmology
  * Delabrouille et al. (2009) Delabrouille, J., Cardoso, J. F., Le Jeune, M., et al. 2009, A&A, 493, 835
  * Dunkley et al. (2013) Dunkley, J. et al. 2013, JCAP, 07, 025
  * Errard et al. (2016) Errard, J., Feeney, S. M., Peiris, H. V., & Jaffe, A. H. 2016, J. Cosmology Astropart. Phys., 2016, 052
  * Errard & Stompor (2012) Errard, J. & Stompor, R. 2012, Phys. Rev. D, 85, 083006
  * Errard & Stompor (2019) Errard, J. & Stompor, R. 2019, Phys. Rev. D, 99, 043529
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759
  * Grain et al. (2009) Grain, J., Tristram, M., & Stompor, R. 2009, Phys. Rev. D, 79, 123515
  * Hamimeche & Lewis (2008) Hamimeche, S. & Lewis, A. 2008, Phys. Rev. D, 77, 103013
  * Harper et al. (2022) Harper, S. E., Dickinson, C., Barr, A., et al. 2022, MNRAS, 513, 5900
  * Harris et al. (2020) Harris, C. R. et al. 2020, Nature, 585, 357
  * Hazumi et al. (2019) Hazumi, M., Ade, P. A. R., Akiba, Y., et al. 2019, Journal of Low Temperature Physics, 194, 443
  * Hervías-Caimapo et al. (2017) Hervías-Caimapo, C., Bonaldi, A., & Brown, M. L. 2017, MNRAS, 468, 4408
  * Hervías-Caimapo et al. (2022) Hervías-Caimapo, C., Bonaldi, A., Brown, M. L., & Huffenberger, K. M. 2022, ApJ, 924, 11
  * Hui et al. (2018) Hui, H., Ade, P. A. R., Ahmed, Z., et al. 2018, Proc. SPIE Int. Soc. Opt. Eng., 10708 [arXiv:1808.00568]
  * Hunter (2007) Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90
  * Ijjas & Steinhardt (2018) Ijjas, A. & Steinhardt, P. J. 2018, Classical and Quantum Gravity, 35, 135004
  * Ijjas & Steinhardt (2019) Ijjas, A. & Steinhardt, P. J. 2019, Physics Letters B, 795, 666
  * Kamionkowski et al. (1997) Kamionkowski, M., Kosowsky, A., & Stebbins, A. 1997, Phys. Rev. Lett., 78, 2058
  * Kamionkowski & Kovetz (2016) Kamionkowski, M. & Kovetz, E. D. 2016, ARA&A, 54, 227
  * Katayama & Komatsu (2011) Katayama, N. & Komatsu, E. 2011, ApJ, 737, 78
  * Krachmalnicoff et al. (2016) Krachmalnicoff, N., Baccigalupi, C., Aumont, J., Bersanelli, M., & Mennella, A. 2016, A&A, 588, A65
  * Krachmalnicoff et al. (2018) Krachmalnicoff, N., Carretti, E., Baccigalupi, C., et al. 2018, A&A, 618, A166
  * Leach et al. (2008) Leach, S. M., Cardoso, J. F., Baccigalupi, C., et al. 2008, A&A, 491, 597
  * Lewis & Challinor (2006) Lewis, A. & Challinor, A. 2006, Phys. Rep, 429, 1
  * Lewis et al. (2000) Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473
  * LiteBIRD Collaboration (2022) LiteBIRD Collaboration. 2022, Progress of Theoretical and Experimental Physics, ptac150
  * Maltoni & Schwetz (2003) Maltoni, M. & Schwetz, T. 2003, Phys. Rev. D, 68, 033020
  * Mangilli et al. (2021) Mangilli, A., Aumont, J., Rotti, A., et al. 2021, A&A, 647, A52
  * Martin et al. (2014a) Martin, J., Ringeval, C., Trotta, R., & Vennin, V. 2014a, J. Cosmology Astropart. Phys., 2014, 039
  * Martin et al. (2014b) Martin, J., Ringeval, C., & Vennin, V. 2014b, Physics of the Dark Universe, 5, 75
  * Namikawa et al. (2022) Namikawa, T., Baleato Lizancos, A., Robertson, N., et al. 2022, Phys. Rev. D, 105, 023511
  * Nash (1984) Nash, S. G. 1984, SIAM Journal on Numerical Analysis, 21, 770
  * Natoli et al. (2018) Natoli, P., Ashdown, M., Banerji, R., et al. 2018, J. Cosmology Astropart. Phys., 2018, 022
  * Planck Collaboration Int. XXX (2016) Planck Collaboration Int. XXX. 2016, A&A, 586, A133
  * Planck Collaboration IV (2020) Planck Collaboration IV. 2020, A&A, 641, A4
  * Planck Collaboration V (2020) Planck Collaboration V. 2020, A&A, 641, A5
  * Planck Collaboration VI (2020) Planck Collaboration VI. 2020, A&A, 641, A6
  * Planck Collaboration X (2016) Planck Collaboration X. 2016, A&A, 594, A10
  * Planck Collaboration X (2020) Planck Collaboration X. 2020, A&A, 641, A10
  * Planck Collaboration XI (2016) Planck Collaboration XI. 2016, Astron. Astrophys., 594, A11
  * Planck Collaboration XI (2020) Planck Collaboration XI. 2020, A&A, 641, A11
  * Poletti & Errard (2023) Poletti, D. & Errard, J. 2023, In prep.
  * Puglisi et al. (2022) Puglisi, G., Mihaylov, G., Panopoulou, G. V., et al. 2022, MNRAS, 511, 2052
  * Remazeilles et al. (2018a) Remazeilles, M., Banday, A. J., Baccigalupi, C., et al. 2018a, J. Cosmology Astropart. Phys., 2018, 023
  * Remazeilles et al. (2011) Remazeilles, M., Delabrouille, J., & Cardoso, J.-F. 2011, MNRAS, 410, 2481
  * Remazeilles et al. (2018b) Remazeilles, M., Dickinson, C., Eriksen, H. K., & Wehus, I. K. 2018b, MNRAS, 474, 3889
  * Remazeilles et al. (2016) Remazeilles, M., Dickinson, C., Eriksen, H. K. K., & Wehus, I. K. 2016, MNRAS, 458, 2032
  * Remazeilles et al. (2021) Remazeilles, M., Rotti, A., & Chluba, J. 2021, MNRAS, 503, 2478
  * Seljak (1997) Seljak, U. 1997, ApJ, 482, 6
  * Seljak & Zaldarriaga (1997) Seljak, U. & Zaldarriaga, M. 1997, Phys. Rev. Lett., 78, 2054
  * Smith & Zaldarriaga (2007) Smith, K. M. & Zaldarriaga, M. 2007, Phys. Rev. D, 76, 043001
  * SO Collaboration (2019) SO Collaboration. 2019, J. Cosmology Astropart. Phys., 2019, 056
  * Starobinskiǐ (1979) Starobinskiǐ, A. A. 1979, Soviet Journal of Experimental and Theoretical Physics Letters, 30, 682
  * Stompor et al. (2016) Stompor, R., Errard, J., & Poletti, D. 2016, Phys. Rev. D, 94, 083526
  * Stompor et al. (2009) Stompor, R., Leach, S., Stivoli, F., & Baccigalupi, C. 2009, MNRAS, 392, 216
  * Tegmark (1998) Tegmark, M. 1998, ApJ, 502, 1
  * Thorne et al. (2019) Thorne, B., Dunkley, J., Alonso, D., et al. 2019, arXiv e-prints, arXiv:1905.08888
  * Thorne et al. (2017) Thorne, B., Dunkley, J., Alonso, D., & Næss, S. 2017, MNRAS, 469, 2821
  * Vacher et al. (2022) Vacher, L., Aumont, J., Montier, L., et al. 2022, A&A, 660, A111
  * Vacher et al. (2023) Vacher, L., Chluba, J., Aumont, J., Rotti, A., & Montier, L. 2023, A&A, 669, A5
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Wagenmakers & Farrell (2004) Wagenmakers, E. & Farrell, S. 2004, Psychon Bull Rev., 11, 192
  * Zaldarriaga & Seljak (1997) Zaldarriaga, M. & Seljak, U. 1997, Phys. Rev. D, 55, 1830
  * Zaldarriaga & Seljak (1998) Zaldarriaga, M. & Seljak, U. 1998, Phys. Rev. D, 58, 023003
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source Software, 4, 1298

##  Appendix A $\hat{\chi}^{2}_{\rm min}$ goodness of fit

Figure 13: Same as Figure 12 but only considering the four less complex
Galactic foreground templates. Note the different $x$ scale. We do not see
strong preference for $C_{\ell}$-moments in any of the cases considered.

Aside from studying the data products at the end of the $B$-mode analysis,
namely the clean CMB power spectra and the marginalized $r$ posterior, we can
take a step back and ask how well the input data are described by our
theoretical models. With pipelines B and C, the input data are multi-frequency
maps that are not directly compared to a likelihood model, instead a map-based
cleaning step produces clean CMB maps that are fitted on the power spectrum
level with a cosmological model. We will therefore not address the goodness of
the model fit with pipelines B and C in this paper. With pipeline A, this
question has a more straightforward answer: here, the input data are multi-
frequency power spectra directly fitted by a $C_{\ell}$-likelihood model
describing CMB, Galactic synchrotron and dust. Another interesting aspect that
can be addressed in the same context is the comparison of the goodness of the
fit between the $C_{\ell}$-fiducial and the $C_{\ell}$-moments model. We will
therefore focus on pipeline A and A + moments. The $\hat{\chi}^{2}_{\rm min}$
statistic is defined as

$\hat{\chi}^{2}_{\rm min}=(\mathbf{d}-\mathbf{t})^{T}{\sf
C}^{-1}(\mathbf{d}-\mathbf{t})\,,$ (17)

where $\mathbf{d}$ is the data vector, $\mathbf{t}$ is the best-fit theory
prediction, and ${\sf C}$ is the data covariance. It can be shown that for
Gaussian distributed data, the minimum-$\chi^{2}$ statistic follows a
$\chi^{2}$ distribution with $N-P$ degrees of freedom, where $N$ is the number
of data points and $P$ is the number of parameters (see e.g. Maltoni & Schwetz
(2003)). In our case ($N-P=567$), this is approximately a Gaussian
distribution with mean $N-P$ and variance $2(N-P)$. If we use Gaussian
foregrounds, our likelihood is accurate since a) the model is known, b) the
power spectrum data are very close to Gaussian (owing to the central limit
theorem, see Section 2.1), and c) the covariance is accurate, since it was
estimated from the corresponding set of simulations. If our guess of the
covariance was an over- or underestimation of the data covariance, then the
$\hat{\chi}^{2}_{\rm min}$ would peak at lower, or higher values,
respectively. We compute the $\hat{\chi}^{2}_{\rm min}$ for 100 simulated data
realizations containing standard cosmology ($r=0$, $A_{\rm lens}=1$),
inhomogeneous noise in the goal-optimistic scenario and in four foreground
cases (Gaussian, d0s0, d1s1, dmsm), considering both models separately. We
observe that the empirical distribution of the $\hat{\chi}^{2}_{\rm min}$
statistic matches the theoretical expectation not only for Gaussian
foregrounds, but also the non-Gaussian cases. This shows that our estimate of
the covariance matrix is appropriate, even for the more complex foreground
models.

While the $\hat{\chi}^{2}_{\rm min}$ statistic evaluates the goodness of the
fit for each model independently, we also make use of a refined model
comparison tool that quantifies the preference toward either model,
considering every individual simulation seed: the Akaike Information Criterion
(AIC, Akaike (1974)). Specifically, we compute the difference $\Delta{\rm
AIC}=2\Delta k+\hat{\chi}^{2}_{\rm min}({\rm moments})-\hat{\chi}^{2}_{\rm
min}({\rm fiducial})$, where $\Delta k=4$ is the number of excess parameters
of the $C_{\ell}$-moments over $C_{\ell}$-fiducial, penalizing the former for
its higher complexity. Following e.g. Wagenmakers & Farrell (2004), the number
$\exp(\Delta{\rm AIC}/2)$ can be interpreted as the relative model odds
$P(C_{\ell}$-fiducial$)/P(C_{\ell}$-moments$)$.

The results are shown in Figure 13. The AIC test detects no clear preference
for $C_{\ell}$-moments over $C_{\ell}$-fiducial. Among 100 simulations
containing dmsm foregrounds, 30 have model odds below 1, and only a single
simulation below $10^{-2}$. Hence, $C_{\ell}$-moments is not clearly preferred
over $C_{\ell}$-fiducial, although it allows us to mitigate a 1-2$\sigma$ bias
on $r$ in the dmsm case. This shows that care must be taken when interpreting
$<2\sigma$ detections in marginal posterior distributions, since they may be
difficult to distinguish from subdominant residual foreground contamination
using standard model-comparison techniques. We will see in Section 5 that this
situation changes in the presence of yet more complex foregrounds that lead to
biases far beyond 2$\sigma$ within the $C_{\ell}$-fiducial model.

##  Appendix B Bias on $r$ for Gaussian foregrounds and pipeline B

Figure 14: Left panel: illustration of foreground residuals in the
reconstructed CMB for inhomogeneous goal-pessimistic noise, with a fiducial
($r=0$, $A_{\rm lens}=1$) CMB, with Gaussian and d0s0 foregrounds. We show the
reconstructed CMB and foreground residual $C_{\ell}^{BB}$. The residual is
calculated by mixing the pure foreground maps with the NILC weights. Right
panel: the more constrained mask (built from the fiducial mask shown in Figure
2) which is closer to homogeneity.

In the results of pipeline B, we note a consistent bias on $r$ for simulations
that include Gaussian foregrounds in combination with inhomogeneous noise.
This bias seems unreasonable, considering the results that use more complex
foregrounds, such as d0s0, but the same noise and CMB simulations, lead to
considerably less bias. The cosmological bias is directly caused by a bias at
the large angular scales of the $BB$ power spectrum, which is easily visible
when plotting the spectra. We confirm that this bias corresponds to foreground
bias. In Figure 14 (left panel), we show the reconstructed CMB $C_{\ell}^{BB}$
spectrum for the standard cosmology ($r=0$, $A_{\rm lens}=1$) with
inhomogeneous goal-pessimistic noise, for both the Gaussian (blue squares) and
d0s0 (red squares) foregrounds. The marker and error bar shows the mean and
1-$\sigma$ standard deviation across the 500 simulations. The excess bias at
large scales is evident as the reconstructed CMB clearly surpasses the
fiducial $BB$ spectrum. We can calculate the exact foreground bias present on
each reconstructed CMB, by mixing the pure foregrounds maps with the same NILC
weights (Equation 5), and taking the power spectra of that foreground bias
reconstruction. This is shown as the blue triangles for Gaussian and red
triangles for d0s0 foregrounds. In the Gaussian case, we see that the large-
scale bias is almost entirely caused by foregrounds.

We know this bias is directly related to the NILC weights that mix the
frequency maps (Equation 5). The weights are calculated directly from the
frequency-frequency covariance matrix, which in turn is calculated from the
input frequency maps themselves. If, for instance, we wish to test if the
weights calculated with the simulations including Gaussian foregrounds are
incorrect, we can do the cross check by using other weights to mix those same
frequency maps. We take the NILC weights calculated for the simulations with
d0s0 foregrounds (shown in Figure 14) and use them to mix the frequency maps
that include the Gaussian foregrounds. In this case, the mean best fit is
$r=-0.0001\pm 0.0026$, which is comparable to $r=-0.0005\pm 0.0028$ for the
simulations including d0s0 foregrounds (shown in Table 4).

Figure 15: Values of the covariance matrix used to build the ILC weights
(Equation 5) from the data. This is plotted for the $B$-field reconstruction,
all frequency-frequency combinations, for a randomly chosen 50% of all pixels.
The left panel shows the first needlet window (largest scales), the right
panel shows the fifth needlet window (smallest scales). Both show the values
for a single simulation seed for inhomogeneous goal-pessimistic noise, with a
fiducial ($r=0$, $A_{\rm lens}=1$) CMB, with Gaussian, d0s0, and d1s1
foregrounds.

Any incorrect weights must originate from the NILC covariance matrix. The
combination of directly calculating the covariance matrix over maps that have
inhomogeneous noise and Gaussian foregrounds creates the observed bias. We can
alleviate this problem by spatially constraining the hits mask (Figure 2) to a
more homogeneous area, where the noise map realization will also be more
homogeneous. The mask we use to do this is shown in Figure 14 (right panel).
This mask has a smaller sky fraction than the fiducial mask we use, and
therefore a price is paid in the form of a greater uncertainty in the
measurement of $r$. The measurements using this more constrained and more
homogeneous mask are marked with a † in Table 4.