of 11
A Simulation-based Method for Correcting Mode Coupling in CMB Angular Power
Spectra
J. S.-Y. Leung
1
,
2
, J. Hartley
3
, J. M. Nagy
4
,
5
, C. B. Netter
fi
eld
1
,
3
, J. A. Shariff
6
, P. A. R. Ade
7
, M. Amiri
8
, S. J. Benton
9
,
A. S. Bergman
9
, R. Bihary
10
, J. J. Bock
11
,
12
, J. R. Bond
6
, J. A. Bonetti
12
, S. A. Bryan
13
, H. C. Chiang
14
,
15
, C. R. Contaldi
16
,
O. Doré
11
,
12
, A. J. Duivenvoorden
9
, H. K. Eriksen
17
, M. Farhang
1
,
6
,
18
, J. P. Filippini
19
,
20
, A. A. Fraisse
9
, K. Freese
21
,
22
,
M. Galloway
17
, A. E. Gambrel
23
, N. N. Gandilo
24
, K. Ganga
25
, R. Gualtieri
26
, J. E. Gudmundsson
22
, M. Halpern
8
,
M. Hassel
fi
eld
27
, G. Hilton
28
, W. Holmes
12
, V. V. Hristov
11
, Z. Huang
6
, K. D. Irwin
29
,
30
, W. C. Jones
9
, A. Karakci
17
,
C. L. Kuo
29
, Z. D. Kermish
9
,S.Li
9
,
31
,D.S.Y.Mak
16
, P. V. Mason
11
, K. Megerian
12
, L. Moncelsi
11
, T. A. Morford
11
,
M. Nolta
6
,R.O
Brient
12
, B. Osherson
19
, I. L. Padilla
1
,
32
, B. Racine
17
, A. S. Rahlin
23
,
33
, C. Reintsema
28
, J. E. Ruhl
10
,
M. C. Runyan
11
, T. M. Ruud
17
, E. C. Shaw
19
, C. Shiu
9
, J. D. Soler
34
, X. Song
9
, A. Trangsrud
11
,
12
, C. Tucker
7
,
R. S. Tucker
11
, A. D. Turner
12
, J. F. van der List
9
, A. C. Weber
12
, I. K. Wehus
17
, S. Wen
10
, D. V. Wiebe
8
, and E. Y. Young
29
,
30
1
Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON, M5S 3H4, Canada;
leung@astro.utoronto.ca
2
Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON, M5S 3H4, Canada
3
Department of Physics, University of Toronto, 60 St. George Street, Toronto, ON, M5S 3H4, Canada
4
Department of Physics, Washington University in St. Louis, 1 Brookings Drive, St. Louis, MO 63130, USA
5
McDonnell Center for the Space Sciences, Washington University in St. Louis, 1 Brookings Drive, St. Louis, MO 63130, USA
6
Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON, M5S 3H8, Canada
7
School of Physics and Astronomy, Cardiff University, The Parade, Cardiff, CF24 3AA, UK
8
Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC, V6T 1Z1, Canada
9
Department of Physics, Princeton University, Jadwin Hall, Princeton, NJ 08544, USA
10
Physics Department, Case Western Reserve University, 10900 Euclid Avenue, Rockefeller Building, Cleveland, OH 44106, USA
11
Division of Physics, Mathematics and Astronomy, California Institute of Technology, MS 367-17, 1200 E. California Boulevard, Pasadena, CA 91125, U
SA
12
Jet Propulsion Laboratory, Pasadena, CA 91109, USA
13
School of Electrical, Computer, and Energy Engineering, Arizona State University, 650 E Tyler Mall, Tempe, AZ 85281, USA
14
Department of Physics, McGill University, 3600 Rue University, Montreal, QC, H3A 2T8, Canada
15
School of Mathematics, Statistics and Computer Science, University of KwaZulu-Natal, Durban, South Africa
16
Blackett Laboratory, Imperial College London, SW7 2AZ, London, UK
17
Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, NO-0315 Oslo, Norway
18
Department of Physics, Shahid Beheshti University, 1983969411, Tehran, Iran
19
Department of Physics, University of Illinois at Urbana-Champaign, 1110 W. Green Street, Urbana, IL 61801, USA
20
Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA
21
Department of Physics, University of Texas, 2515 Speedway, C1600, Austin, TX 78712, USA
22
The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, SE-106 91 Stockholm, Sweden
23
Kavli Institute for Cosmological Physics, University of Chicago, 5640 S Ellis Avenue, Chicago, IL 60637, USA
24
Steward Observatory, 933 North Cherry Avenue, Tucson, AZ 85721, USA
25
APC, Univ. Paris Diderot, CNRS
/
IN2P3, CEA
/
Irfu, Obs de Paris, Sorbonne Paris Cité, France
26
High Energy Physics Division, Argonne National Laboratory, Argonne, IL 60439, USA
27
Department of Astronomy and Astrophysics, Pennsylvania State University, 520 Davey Lab, University Park, PA 16802, USA
28
National Institute of Standards and Technology, 325 Broadway Mailcode 817.03, Boulder, CO 80305, USA
29
Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA
30
SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA
31
Department of Mechanical and Aerospace Engineering, Princeton University, Engineering Quadrangle, Princeton, NJ 08544, USA
32
Department of Physics and Astronomy, Johns Hopkins University, 3701 San Martin Drive, Baltimore, MD 21218, USA
33
Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510-5011, USA
34
Max-Planck-Institute for Astronomy, Konigstuhl 17, D-69117, Heidelberg, Germany
Received 2021 November 2; revised 2022 February 1; accepted 2022 February 11; published 2022 March 30
Abstract
Modern cosmic microwave background
(
CMB
)
analysis pipelines regularly employ complex time-domain
fi
lters,
beam models, masking, and other techniques during the production of sky maps and their corresponding angular
power spectra. However, these processes can generate couplings between multipoles from the same spectrum and
from different spectra, in addition to the typical power attenuation. Within the context of pseudo-
C
based,
MASTER
-style analyses, the net effect of the time-domain
fi
ltering is commonly approximated by a multiplicative
transfer function,
F
, that can fail to capture mode mixing and is dependent on the spectrum of the signal. To
address these shortcomings, we have developed a simulation-based spectral correction approach that constructs a
two-dimensional transfer matrix,
J
ℓℓ
¢
, which contains information about mode mixing in addition to mode
attenuation. We demonstrate the application of this approach on data from the
fi
rst
fl
ight of the
SPIDER
balloon-
borne CMB experiment.
The Astrophysical Journal,
928:109
(
11pp
)
, 2022 April 1
https:
//
doi.org
/
10.3847
/
1538-4357
/
ac562f
© 2022. The Author
(
s
)
. Published by the American Astronomical Society.
Original content from this work may be used under the terms
of the
Creative Commons Attribution 4.0 licence
. Any further
distribution of this work must maintain attribution to the author
(
s
)
and the title
of the work, journal citation and DOI.
1
Uni
fi
ed Astronomy Thesaurus concepts:
Cosmic microwave background radiation
(
322
)
;
Observational cosmology
(
1146
)
;
Astronomy data analysis
(
1858
)
1. Introduction
Producing well-characterized maps of the cosmic microwave
background
(
CMB
)
from time-ordered data requires accurately
accounting for the impact of instrumental effects and any signal
processing on the underlying astrophysical signal
(
e.g., Jarosik
et al.
2007
; Planck Collaboration et al.
2016
)
. These processing
techniques typically also have a nontrivial impact on the
fi
delity of the cosmological signal in ways that are spatially
anisotropic and inhomogeneous. These effects must be
precisely and accurately characterized in order to avoid biasing
the estimation of the CMB angular power spectra, and therefore
the inference of cosmological parameters.
The impact of timestream processing and instrument
response is commonly approximated in the harmonic domain
with a
fi
lter window function unique to the instrument, derived
from the analysis of large ensembles of signal simulations. In
its simplest formulation, the processing pipeline is modeled as
a power attenuation mechanism in each multipole, i.e., the
fi
lter
window function is a transfer function
a simple one-
dimensional vector of ratios of output to input power.
This paper addresses the shortcomings of the one-dimen-
sional model and proposes an alternative approach through the
construction of a two-dimensional transfer matrix. This
simulation-based spectral correction approach takes into
account both the mode mixing and attenuation from instru-
mental and data processing effects including beams,
fi
ltering,
and masking. Section
2
of this paper describes the theoretical
motivation and compares the transfer matrix to the one-
dimensional transfer function approach. A concrete example is
presented in Section
3
using data from the
fi
rst
fl
ight of
SPIDER
,
a balloon-borne telescope designed to measure CMB polariza-
tion on roughly degree angular scales
(
SPIDER Collabora-
tion
2021
)
. This section explores different techniques for
constructing the transfer matrix to reduce the computational
demand including using binned power spectra and performing
Fourier-space interpolation. Section
4
presents several compar-
isons of these techniques, including tests of signal recovery on
spectra with different shapes. While all tested approaches were
found to accurately recover a target spectrum identical to that
used for the transfer matrix construction, the performance
varied when applied to a different target spectrum. As
discussed in Section
5
, this has implications for increasingly
sensitive CMB polarization measurements where the cosmolo-
gical signals are heavily obscured by Galactic foregrounds. As
the foreground power spectra are less well constrained and vary
substantially between different sky regions, understanding the
signal dependence of potential analysis techniques becomes
even more important.
2. Theoretical Description
Maps of the CMB temperature
(
T
)
and linear polarization
(
Q
,
U
)
anisotropies over the partial sky can be decomposed into
linear combinations of spherical harmonics:
rrr
TW
aY
,1
ℓm
ℓm
T
ℓm
,
() ()
̃
()
( )
å
=
rrrr
QiUW aY
,2
ℓm
ℓm
ℓm
,
22
[()
()] ()
̃
()
()
å
=

where
W
(
r
)
is the window representing the relative weights of
the partial-sky mask. To avoid using spin-weighted spherical
harmonic components, the
a
ℓm
2
̃
coef
fi
cients are frequently
expressed as a combination of scalar and
fi
xed parity
E
and
B
components:
aaia
,3
ℓm
ℓm
E
ℓm
B
2
̃( ̃ ̃)
()
=- 
where the sign convention follows Zaldarriaga & Seljak
(
1997
)
. For each
a
ℓm
̃
, the pseudo-power spectrum
C
̃
is de
fi
ned
as
C
a
1
21
.4
m
ℓm
2
̃
∣ ̃ ∣
( )
å
=
+
Also known as a pseudo-
C
(
PCL
)
, it is related to the angular
power spectrum speci
fi
ed by the theory of primordial
perturbations,
C
, via
CKC
,5
X
ℓX
ℓℓ
XX
X
,
̃
()
å
áñ=
¢¢
¢
¢
¢
¢
where
K
ℓℓ
XX
¢
¢
is a mode coupling kernel
(
or
mixing matrix
)
that
accounts for the mixing within and between the observables
X
ä
{
TT
,
EE
,
BB
,
TE
,
EB
,
TB
}
due to the partial-sky mask, and
the brackets
·
denote an ensemble average over in
fi
nite
spectrum realizations. Because
K
ℓℓ
XX
¢
¢
is entirely determined by
the chosen pixel weighting and the geometry of the cut sky
(
Hivon et al.
2002
)
, in the absence of instrumental effects and
noise, PCL estimators can use the
(
fi
nite
)
set of measured
C
̃
to
solve for the underlying power spectrum
C
. Examples of PCL
estimators include
MASTER
(
Hivon et al.
2002
)
,
NaMaster
(
Alonso Monge et al.
2019
)
, and
PolSpice
(
Chon et al.
2004
)
.
A major challenge in interpreting CMB data is that the
instruments cannot directly probe the true sky signal; the
incoming signals are inevitably altered by instrumental
systematics and noise. Therefore, as described in Section
1
,
an experiment
s raw datastreams must be processed to remove
many different types of spurious signals. The act of observing
and time-domain
fi
ltering both distort the signal estimate and
present additional sources of mode coupling. Assuming that the
coupling is homogeneous, isotropic, and linear, the impact of
the experiment is captured by introducing another coupling
matrix,
F
ℓℓ
XX
¢
¢
, such that
CKFC
.6
X
ℓXℓX
ℓℓ
XX
ℓℓ
XX
X
,,
̃
()
åå
áñ=
¢¢
¢
¢
¢
¢
In general,
F
ℓℓ
XX
¢
¢
is unique to the experiment and cannot be
determined analytically.
2.1. The Filter Transfer Function
To reduce complexity,
F
ℓℓ
XX
¢
¢
is often approximated as a
diagonal matrix whose entries represent the one-dimensional
2
The Astrophysical Journal,
928:109
(
11pp
)
, 2022 April 1
Leung et al.
transfer function
F
XX
¢
, such that
CKFC
7
ℓℓ ℓ ℓ
̃
()
å
áñ=
¢
¢¢¢
(
here and hereafter, the superscripts
X
are suppressed for
brevity
)
. Using a set of pseudo-power spectra
C
̃
derived from
an ensemble of simulations of a known power spectrum
C
, the
transfer function
F
can be estimated through an iterative
process to avoid inverting
K
ℓℓ
¢
(
Hivon et al.
2002
; Dutcher et al.
2021
)
. Note that
F
is frequently decomposed further into a
fi
lter component
f
and a beam component
b
2
, such that
Ffb
;
2
=
additional instrumental effects can be inserted in a
similar fashion.
This one-dimensional approximation implicitly assumes that
each
-mode remains independent throughout the entire
fi
ltering process. As long as this assumption holds, the
mapping from input to output is one-to-one:
F
is simply the
ratio of output to input power for each
. However, in practice,
we expect modes to become coupled with one another, where
the mapping becomes many-to-one and the contributions from
the coupled input modes become impossible to disentangle
using
F
alone. The consequence is that changing the input
power spectrum
C
also changes the output
C
̃
in some
nontrivial way due to this many-to-one mapping. In other
words,
F
is inextricably tied to the input used to compute it.
More concretely, the one-dimensional transfer function
formulation con
fl
ates mode mixing with the in-mode
fi
lter
gain. This introduces a sensitivity to the power spectrum of the
simulated sky used to calibrate
F
; the
fi
nal spectrum is correct
only if the simulated sky is statistically similar to the true sky
(
i.e., a Gaussian sky realization with the same power spectra
)
.
2.2. The Multipole
Multipole Transfer Matrix
To address this shortcoming of the
F
formulation, we
introduce a two-dimensional linear operator that encodes the
(
asymmetric
)
coupling between each
ℓℓ
¢
pair:
CJC
.8
ℓℓ ℓ
̃
()
å
áñ=
¢
¢¢
We refer to this coupling operator
J
ℓℓ
¢
as the
transfer matrix
because it directly relates the input of the true power spectrum
to the output of the spectrum estimator. This relation holds as
long as the
fi
ltering process is approximately linear. Treating
the entire pipeline as a single operator avoids the diagonal
approximation and ensures that all multipole
multipole
couplings induced by time-domain and map-domain
fi
ltering
are properly taken into account, i.e., we are not locked into a
speci
fi
c input spectrum.
Note that these couplings extend to those between the six
power spectra; both temperature-to-polarization
(
T
-to-
P
)
and
E
-
to-
B
leakage are automatically included. Because standard PCL
spectrum estimators provide the
TT
,
EE
,
BB
,
TE
,
EB
, and
TB
power spectra, any input mode can be readily related to an
output mode from any of these six power spectra, resulting in
6
2
coupling matrices. We
fi
nd it convenient to compile these
individual matrices into a single 6
×
6 block matrix encapsulat-
ing every
ℓℓ
¢
coupling between the six power spectra.
Following the rules of matrix multiplication
(
Equation
(
8
))
,
we arrange these blocks horizontally according to input spectra
and vertically according to output spectra
(
see Figure
1
)
. The
ordering of the six spectra does not matter as long as it is
consistent between the two axes; likewise, the six
C
s must be
concatenated in the same order. We choose the above ordering
based on convenience.
While the approach presented here has similarities to that
used by the BICEP
/
Keck Collaboration for their CMB
polarization analyses, the implementation varies due to key
differences in the observing strategies and analysis pipelines.
As described in Ade et al.
(
2016
)
, the simplicity of the BICEP
/
Keck observing strategy allows for the construction of an
observing matrix that renders large simulation ensembles more
computationally tractable. Having determined their
fi
ltering
operations to be linear, the observing matrix is used to
construct the transfer matrix
J
ℓℓ
¢
, which is then used similarly to
reconstruct and interpret the on-sky power spectra. Because the
observing matrix formulation is less generally applicable to
experiments at other sites with different observing strategies,
the remainder of this paper is dedicated to estimating
J
ℓℓ
¢
through other means.
3. Application to
SPIDER
Data
To illustrate the application of the transfer matrix, we present
results from the
SPIDER
balloon-borne telescope. During its
fi
rst Antarctic long-duration balloon
(
LDB
)
fl
ight in January
2015,
SPIDER
mapped 4.8% of the sky with polarimeters
operating at 95 and 150 GHz to constrain the
B
-mode power
spectrum from primordial gravitational waves
(
SPIDER
Collaboration
2021
)
. An upcoming
fl
ight with additional
280 GHz receivers will provide improved characterization of
the polarized Galactic dust foregrounds
(
Shaw et al.
2020
)
.
SPIDER
s processing pipeline is described more fully in
SPIDER Collaboration
(
2021
)
and SPIDER Collaboration
(
2022, in preparation
)
; here we brie
fl
y highlight the most
relevant steps and note that they are suf
fi
ciently linear to allow
usage of Equation
(
8
)
. Once features such as cosmic-ray hits,
payload transmitter signals, and thermal transients have been
removed from the raw detector timestreams, they are
fi
ltered to
reduce quasi-stationary noise. Null test performance was used
to identify the weakest
fi
lter that suf
fi
ciently removed quasi-
stationary noise, resulting in a
fi
fth-order polynomial
fi
t to each
detector
s data as a function of azimuth angle between scan
turnarounds. The impact of scanning,
fi
ltering, and
fl
agging are
determined by applying the entire processing pipeline to an
ensemble of time-domain signal simulations in a procedure
known as reobservation. The reobserved timestreams are
produced at the full data sample rate without any down-
sampling or binning. Unlike the BICEP
/
Keck experiment,
SPIDER
s observing strategy does not allow for the creation of
an observing matrix, so obtaining the transfer matrix
J
ℓℓ
¢
requires producing an appropriate set of reobserved
CMB maps.
SPIDER
s measured and simulated timestreams are converted
into two-dimensional maps of the sky with a binned mapmaker.
This approach assembles the detector data into spatial pixels
based on the telescope pointing and polarization sensitivity as
described further in SPIDER Collaboration
(
2021
)
. The
computational simplicity of this method is crucial for enabling
the generation of large simulation ensembles, including those
used in this work.
The main
SPIDER
cosmological results presented in SPIDER
Collaboration
(
2021
)
use two complementary pipelines for
power spectrum estimation. Here we describe results only from
the Noise Simulation Independent
(
NSI
)
pipeline, while the
3
The Astrophysical Journal,
928:109
(
11pp
)
, 2022 April 1
Leung et al.
other pipeline is presented in Gambrel et al.
(
2021
)
. The NSI
pipeline is similar to Xspect
(
Tristram et al.
2005
)
and Xpol
(
Tristram
2006
)
, as well as to those used by the South Pole
Telescope
(
Lueker et al.
2010
)
, CLASS
(
Padilla et al.
2020
)
,
and
SPIDER
s circular polarization analysis
(
Nagy et al.
2017
)
.
It decomposes the 2015
fl
ight data into 14 maps by splitting the
timestreams into interleaved 3 minute segments for each of the
six receivers. This timescale was chosen to maximize the
number of complete and well-conditioned maps of the sky
region, taking into account
SPIDER
s scan rate and observing
strategy. These maps are then co-added by frequency to
produce 14 maps in each of the two
SPIDER
bands, 95 and
150 GHz. The cross-power spectra for all pairs of maps
(
neglecting the auto-spectra
)
are produced with
PolSpice
(
Chon et al.
2004
)
, and the distribution of the cross-spectra
allows for an empirical measurement of the uncertainty without
reliance on an instrumental noise model. In
SPIDER
s NSI
pipeline, these cross-spectra are used to estimate cosmological
parameters by comparing them to theoretical models that
depend on the parameters of interest. It is therefore
SPIDER
s
power spectra, rather than the maps, that we correct for the
beam and
fi
ltering effects as described in the rest of this work.
Since the
PolSpice
estimator is linear, the assumption of
linearity in the de
fi
nition of the transfer matrix
(
Equation
(
8
))
is
satis
fi
ed. Although the transfer matrix can be used to correct
the power spectra for temperature-to-polarization leakage, we
perform this correction on the maps based on simulations of
Planck temperature-only maps
(
Planck Collaboration et al.
2020
)
. This avoids the sample variance from the simulation
ensemble, which is relatively large due to the amplitude of the
temperature signals, by using existing measurements of the
temperature anisotropies in
SPIDER
s observing region.
As
J
ℓℓ
¢
is merely the linear operator linking the input and
output power spectra
(
Equation
(
8
))
, its computation is
conceptually straightforward: pass a
(
simulated
)
CMB map
with a known spectrum through the reobservation pipeline,
then compute the power spectra of the output map. But to
capture the asymmetric mode
mode coupling between each
ℓℓ
¢
pair, this process must be performed on each mode individu-
ally. Thus we begin with a set of unit
δ
-functions
one at each
of interest
and simulate a set of CMB maps using the
synfast HEALPix
utility
35
(
Gorski et al.
2005
)
. Upon
obtaining the reobserved maps, we organize their power spectra
into columns, arranged by
, to form
J
ℓℓ
¢
. The whole process is
illustrated schematically in Figure
1
.
To reduce the effects of sample variance from the CMB map
realization and partial-sky spectrum estimation procedures, we
repeat the above process 100 times and average the 100
resulting matrices to obtain the
fi
nal
J
ℓℓ
¢
. However, the
computation cost is intensive: with 100 map realizations of a
δ
-function in each of the three CMB modes
(
TT
,
EE
,
BB
)
that
must be reobserved by each of
SPIDER
s six receivers, we
require 1800 simulated maps per
ä
[
2, 250
+
Δ
b
]
(
a
Δ
b
buffer
is necessary to account for couplings of higher
multipoles with lower multipoles
)
. To minimize the computa-
tional burden, we use full-mission
SPIDER
maps instead of the
14 NSI pipeline submaps because this was demonstrated to
have a negligible effect on
SPIDER
s signal recovery in
simulations. Nevertheless, reobserving the entire set of required
maps would take about 1000 core-years on the
Niagara
supercomputing cluster
(
Loken et al.
2010
; Ponce et al.
2019
)
,
which is computationally infeasible given
SPIDER
Collabora-
tion resources. Since the
fi
delity of the computation needs to be
Figure 1.
Procedure to create the transfer matrix
J
ℓℓ
¢
. For each input mode, the six output response spectra
(
gains
)
are arranged vertically to form the columns of the
matrix. Each of the 36 blocks within the matrix spans a range of multipoles with

ℓℓℓ
min
max
.
35
https:
//
healpix.sourceforge.io
4
The Astrophysical Journal,
928:109
(
11pp
)
, 2022 April 1
Leung et al.
balanced against the time required, two practical options are
investigated in the following sections:
(
1
)
use input spectra
with multipole bins rather than evaluating each multipole
individually, or
(
2
)
compute the transfer matrix for only a
sample of input
s and interpolate between them.
3.1. Bin
Bin Transfer Matrix
CMB angular power spectra are typically binned into
multipole ranges of width
Δ
in order to increase the signal-
to-noise ratio and mitigate the correlative effects of adjacent
multipole leakage. As a result, we consider the case where the
input power spectra are also binned into the same multipole
ranges, thereby reducing the computational demands of the
calculation. By analogy with Equation
(
8
)
,wede
fi
ne a bin
bin
transfer matrix
J
bb
¢
that describes the average power in output
bin
b
given the average power in each of the input bins
b
¢
:
CJPCJC
,9
b
b
bb
bℓ ℓ
b
bb b
̃
()
åå å
áñ=
=
¢
¢¢
¢
¢¢
where
P
b
is a binning operator.
At
fi
rst glance, the
δ
-function inputs to
J
ℓℓ
¢
might be replaced
by unit-boxcars for
J
bb
¢
, but, as it turns out, that may not be the
best choice. Unlike the case for the full
J
ℓℓ
¢
, the binning
introduces a dependence on the shape of the input spectrum
used to estimate the transfer matrix
(
see the
Appendix
)
.
Consequently, the choice of input spectrum should be as close
as possible to the anticipated output spectrum, as is further
discussed in Section
4
. Because our application is toward
degree-scale CMB
B
-modes, we use a
Λ
cold dark matter
(
CDM
)
+
(
r
=
0.03
)
spectrum, provided by
CAMB
(
Lewis et al.
2000
)
, as our base model. Windowing this model spectrum at
each bin of interest
i.e., multiplying it by unit-boxcars
centered on each bin
provides the set of input spectra needed
to construct the bin
bin transfer matrix
J
bb
¢
.
Because the analysis of
SPIDER
data uses only 12 bins in the
range 2
„
„
300, the computational demand is decreased by a
factor of
30; months of cluster time is reduced to days. The
price is a dependence on an input model, but unlike
F
, the
ability of
J
bb
¢
to capture couplings between multipole bins
provides a more general framework.
3.2. Fourier-space Interpolation
An alternative method for reducing computation time is to
compute the transfer matrix
J
ℓℓ
¢
only at select multipoles
(
i.e.,
columns
)
and interpolate in between. This requires the matrix
to be smooth in the
-range of interest, but, as we shall see, this
assumption is often well justi
fi
ed in practice.
Because
J
ℓℓ
¢
is a linear operator
(
Equation
(
8
))
, the output
from a
δ
-function input can be interpreted as the impulse
response of a linear
fi
lter. As shown in the left panel of
Figure
2
,
SPIDER
s transfer components
(
TT
TT
,
EE
EE
,
etc.
)
have impulse responses resembling displaced sinc
functions. This allows us to recast the problem from the
domain to its dual frequency-like domain
(
hereafter simply the
frequency
domain
)
, where the interpolation turns out to be
considerably easier. For our application, we use the sign and
normalization conventions established by NumPy
s
fft
module;
36
for a given input
¢
, the discrete Fourier transform
of its impulse response
x
the frequency response
is de
fi
ned
by
Xxi
f
N
exp
2
.10
f
N
0
1
()
å
p
=-
=
-
We use
N
=
801 points
(
0
„
„
800
)
in our computations. The
magnitudes
|
X
f
|
and phases
f
=
X
f
of the resulting transforms
are shown in the right panel of Figure
2
.
In general, the properties of the responses are unique to each
experiment and therefore often determined empirically.
SPI-
DER
s response corresponds to a linear-phase low-pass
fi
lter:
the magnitude rolls off until it hits the noise
fl
oor of the
fi
lter at
f
=
f
c
, while the phase
f
is nearly linear up to that same
frequency. The slope of the linear-phase component is
proportional to the shift
(
in
)
of the impulse, namely,
m
=
f
/
f
=
2
π
when
f
is expressed in normalized frequency
units
(
cycles per multipole
)
. The magnitude of the noise
fl
oor
encodes the level of background noise at the input multipole
¢
.
This simple structure in the frequency domain permits a
straightforward interpolation. The magnitudes, in particular,
can be interpolated directly at all
f
, as can the unwrapped
phases at
f
<
f
c
. For simplicity, we interpolate piecewise-
linearly. In the noise region
f
>
f
c
, the
wrapped
phases exhibit a
curious behavior: they settle onto one of two possible values
separated by
π
. Thus, rather than interpolating, we implement a
nearest-neighbor scheme for the phases at
f
>
f
c
.
One might be tempted to
fi
lter out the noise at
f
>
f
c
completely. This requires care: any operation in the frequency
domain carries consequences into the
domain. For example, a
simple low-pass
fi
lter with a rigid cutoff is equivalent to
convolution with a sinc function in
-space; this operation will
pollute the impulse response with heavy ringing, and should be
avoided. Because good
fi
lter design merits a study of its own,
we do not carry out any
fi
ltering of this sort.
To complete the process, the interpolated values are inverse
Fourier transformed back into
-space to form our interpolated
J
ℓℓ
interp
¢
(
Figure
3
)
. The choice of nodes for the interpolation, i.e.,
which columns of
J
ℓℓ
¢
to compute rather than interpolate,
reduces to a balance between
fi
delity and computation time.
Naturally we will want to concentrate our efforts to regions
where
J
ℓℓ
¢
changes rapidly in order to obtain a good
interpolation. For us, that is the low-
regime. We thus
compute the columns at
=
{
5, 10, 15
}
in addition to
SPIDER
s
bandcenters.
We note that this particular choice of interpolation method is
motivated by
SPIDER
s sinc-function-like impulse responses.
However, any suf
fi
ciently smooth
J
ℓℓ
¢
can be interpolated with
similar techniques.
3.3. Choice of Spectrum Estimator
The matrix-making process encodes the spectrum estimator
into the matrix
that is,
J
ℓℓ
¢
depends on the choice of
estimator by design
(
Equation
(
8
))
. Nonetheless, this
choice can be trivially changed because the spectrum
estimation and reobservation
steps are completely indepen-
dent of each other.
The
anafast
utility
(
provided in the
healpy
package
)
is
a simple spectrum estimator that directly computes the power
spectra from a set of input maps using Equation
(
4
)
. In general,
the interpolated matrix constructed using
anafast
provides
very good agreement with the expected results. However, we
fi
nd that the
buffer
needs to be quite large
(
Δ
b
>
120
)
in
36
https:
//
numpy.org
/
doc
/
stable
/
reference
/
routines.fft.html
5
The Astrophysical Journal,
928:109
(
11pp
)
, 2022 April 1
Leung et al.
order to properly account for
E
-to-
B
leakage
(
Figure
4
)
.To
avoid the additional computational demand, in the following
work, we proceed using the
PolSpice
estimator, which
includes as part of its algorithm a correction for
E
-to-
B
leakage
built-in. This reduces the
buffer
to a much more manageable
Δ
b
=
50.
We use the same
PolSpice
parameters as the
SPIDER
B
-
mode analysis for this work, including a cosine apodization
with
σ
=
10
°
,
50
max
q
=
, and
symmetric
_
cl
(
SPIDER
Collaboration
2021
)
.
3.4. Cross-spectra: TE, EB, and TB
An outstanding problem in our discussion so far is the
treatment of the cross-spectra
TE
,
EB
, and
TB
. While the inter-
spectral couplings from auto- to cross-spectra
(
TT
TE
, etc.
)
Figure 2.
Left: impulse responses of unit
δ
-function inputs at
=
20, 45,...,295 for the 150 GHz
EE
EE
coupling; the
TT
and
BB
counterparts are similar. Right: the
corresponding magnitude and
(
unwrapped
)
phase responses. The
frequency
axis has units cycles
/
and does not have a simple physical interpretation. These
responses are interpolated in the Fourier domain to approximate the transfer matrix
(
Figure
3
)
.
Figure 3.
The impulse responses from the left panel of Figure
2
, arranged into columns
(
left
)
, are interpolated in Fourier space to form an approximated
J
ℓℓ
interp
¢
(
right
)
.
Shown is the 150 GHz
EE
EE
coupling; the
TT
and
BB
counterparts are similar.
6
The Astrophysical Journal,
928:109
(
11pp
)
, 2022 April 1
Leung et al.
are implicitly accounted for by the spectrum estimator
(
which
returns all six power spectra
)
, the method described here cannot
be used to compute the matrix components with a cross-
spectrum as input. To wit, it is impossible to generate maps of
Gaussian-distributed
a
m
s with a nonzero correlation
(
the input
mode
)
in
TE
,
EB
,or
TB
while the
TT
,
EE
, and
BB
spectra are
imposed to be zero. Consequently, we approximate the
diagonal cross-spectra transfer components
(
TE
TE
, etc.
)
as
the geometric mean of the related auto-spectra components,
JJJ
ℓℓ
TE TE
ℓℓ
TT TT
ℓℓ
EE EE
,,,
12
(∣
∣∣
∣)
=
¢¢¢
, etc., and treat the off-diagonal
components
(
TE
EB
, etc.
)
as negligible. Although this is not
ideal, it was found not to bias simulation results for our
application. An alternative method could be to construct the
transfer matrix using the harmonic modes, rather than the
power spectra, of the maps. As an example, see Choi et al.
(
2020
)
, who postulate a transfer matrix built from
(
two-
dimensional
)
Fourier transforms of
T
,
E
, and
B
maps.
4. Comparison of Methods
We quantify the effectiveness of the two methods described
above
(
binning and interpolating
)
by two metrics: their ability
to replicate the end-to-end
SPIDER
pipeline, and to recover the
original signal. In the context of Equation
(
8
)
, this corresponds
to producing the correct
C
̃
when given
C
, and vice versa.
4.1. Case 1: Solve for
C
̃
, Given
C
We demonstrate the
fi
rst of these by applying our computed
transfer function
/
matrices directly to a theory spectrum. This
result is then compared to the output spectrum generated by
passing a set of map realizations of the same input theory
spectrum through the
SPIDER
reobservation procedure. To
ensure broad applicability, we use two qualitatively different
theory spectra as input: a
Λ
CDM
+
(
r
=
0.03
)
spectrum
(
from
CAMB
; 500 realizations
)
and a power-law dust foreground
spectrum
(
300 realizations
)
. These two spectra are chosen for
their contrasting behavior
(
Figure
5
)
.
As discussed in Section
2.1
, the input-dependence of
F
makes it suitable only when the target spectrum does not
deviate far from the input spectrum that was used to construct
F
. This is demonstrated in Figure
6
. Having derived our
F
using a
fi
ducial
Λ
CDM
+
(
r
=
0.03
)
spectrum as input,
applying it to the same spectrum again outputs a trivial result
(
left panels
)
. But when the target is the foreground spectrum
(
right panels
)
, large deviations occur wherever mode
mode
coupling dominates over in-mode gain
i.e., wherever the
transfer matrix is least diagonal-like. In our case, this occurs at
low
s
(
Figure
3
)
. The simple power attenuation model of
F
is
insuf
fi
cient in capturing
SPIDER
s
fi
ltering scheme in this
regime.
The same kind of input-dependence is seen in the binned
transfer matrix
J
bb
¢
(
see the
Appendix
)
. Like
F
, the
J
bb
¢
constructed from a
fi
ducial
Λ
CDM
+
(
r
=
0.03
)
spectrum is
highly effective when the target spectrum is also a
Λ
CDM
+
(
r
=
0.03
)
spectrum, and diverges when the target
is the foreground dust spectrum
(
Figure
6
)
. But this time, all
foreground bins up to
=
100 are noticeably af
fl
icted. This is
due to the loss of granularity from binning the spectrum and the
matrix, as we had used a bin width of
Δ
=
25 to match
SPIDER
s
B
-mode analysis. In the limit where the bin width is
1
a single multipole
we recover the full, theoretically ideal
transfer matrix
J
ℓℓ
¢
.
The interpolated transfer matrix
J
ℓℓ
interp
¢
thus functions as a
compromise allowing us to have a bin width of 1 without being
computationally intractable. Despite the simple piecewise-
linear interpolation approximation used, it reliably recovers
both target spectra.
4.2. Case 2: Solve for
C
, Given
C
̃
To recover
C
from
C
̃
, the linear system de
fi
ned by
Equation
(
8
)
must be solved numerically; directly inverting
Figure 4.
Although
anafast
does not automatically correct for
E
-to-
B
leakage, a transfer matrix built using this utility can still capture such mode
mixing induced by
SPIDER
s
fi
lter and mask. When the matrix is applied to a
Λ
CDM
+
(
r
=
0.03
)
theory spectrum, the result can predict the output from
reobserving a map with the same theory spectrum
(
blue points; error bars
indicate the error on the mean of the 500 simulations
)
. However, replicating the
full effects of the reobservation procedure in the higher bins requires
computing a larger matrix than anticipated due to the wide tails of the
EE
BB
response. Namely, a matrix with
300
max
=
(
dotted line
)
under-
estimates the power in the higher bins; this is alleviated somewhat by
increasing to
370
max
=
(
solid line
)
, but that is still not quite enough. As a
result, we proceed with a matrix constructed using the
PolSpice
estimator, as
its built-in
E
-to-
B
leakage correction algorithm avoids the need to generate
larger matrices, thereby reducing the computational demand.
Figure 5.
Two
fi
ducial input theory spectra with very different shapes are used
to evaluate the performance of different methods of computing the transfer
function: a theoretical
Λ
CDM CMB spectrum and a power-law Galactic
foreground model. Here only the CMB
EE
spectrum is shown, while the same
foreground model is used for temperature and polarization.
7
The Astrophysical Journal,
928:109
(
11pp
)
, 2022 April 1
Leung et al.
a linear system is not recommended except in special cases.
One such case is the transfer function
F
computing its
inverse is trivial since it is a diagonal matrix
(
a one-
dimensional function
)
. Likewise,
SPIDER
s binned transfer
matrix
J
bb
¢
is also easily invertible because it is diagonally
dominant and measures only 72
×
72
(
12 bins in each of the
six correlation spectra
)
; its condition number is very close to
1. The main challenge, therefore, is solving the system for the
interpolated transfer matrix
J
ℓℓ
interp
¢
, whose condition number
is
10
17
()
.
The large condition number of
J
ℓℓ
interp
¢
and also of
J
ℓℓ
¢
in
general
means that it is highly susceptible to numerical noise
Figure 6.
Average ratios of reobserved 150 GHz spectra to those obtained from applying the transfer function and transfer matrices to
fi
ducial
Λ
CDM
+
(
r
=
0.03
)
and foreground input spectra
(
500 and 300 realizations, respectively
)
. A ratio of 1 represents the ideal performance, signifying a perfect replication of
SPIDER
s
fi
ltering process. As the transfer function
F
was derived using the reobserved
Λ
CDM result, applying it to that same spectrum again returns 1 by de
fi
nition
(
left
panels
)
. Error bars indicate the error on the mean of the simulation ensembles.
Figure 7.
Solutions to the linear system
Jx JC
ℓℓ
ℓℓ
interp
interp theory
å
=
å
¢
¢
¢
¢
¢
for a
Λ
CDM
EE
spectrum
(
left
)
and a dust foreground spectrum
(
right
)
computed by the
truncated singular value decomposition
(
TSVD
)
method with 50 singular values and the biconjugate gradient stabilized
(
BiCGSTAB
)
method. Near-optimal results
can be obtained when using the TSVD matrix as a preconditioner and providing an initial guess close to
C
theory
(
details provided in the text
)
. Though not shown here,
TT
and
BB
spectra return similar results.
8
The Astrophysical Journal,
928:109
(
11pp
)
, 2022 April 1
Leung et al.