COMAP Pathfinder – Season 2 results - I. Improved data selection and processing

A&A, 691, A335 (2024)

https:

doi.org

10.1051

0004-6361

202451121

The Authors 2024

Astronomy

Astrophysics

COMAP Pathfinder – Season 2 results

I. Improved data selection and processing

J. G. S. Lunde

, N.-O. Stutzer

, P. C. Breysse

6,7

, D. T. Chung

3,4,5

, K. A. Cleary

, D. A. Dunne

H. K. Eriksen

, S. E. Harper

, H. T. Ihle

, J. W. Lamb

, T. J. Pearson

, L. Philip

12,16

, I. K. Wehus

D. P. Woody

, J. R. Bond

3,14,15

, S. E. Church

, T. Gaier

, J. O. Gundersen

, A. I. Harris

, R. Hobbs

J. Kim

, C. R. Lawrence

, N. Murray

3,14,15

, H. Padmanabhan

, A. C. S. Readhead

T. J. Rennie

10,11

, and D. Tolgay

3,14

(COMAP Collaboration)

Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029, Blindern N-0315, Oslo, Norway

California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA

Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada

Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto ON M5S 3H4, Canada

Department of Astronomy, Cornell University, Ithaca, NY 14853, USA

Center for Cosmology and Particle Physics, Department of Physics, New York University, 726 Broadway, New York, NY 10003,

USA

Department of Physics, Southern Methodist University, Dallas TX 75275, USA

Departement de Physique Théorique, Universite de Genève, 24 Quai Ernest-Ansermet, CH-1211 Genève 4, Switzerland

Owens Valley Radio Observatory, California Institute of Technology, Big Pine, CA 93513, USA

Jodrell Bank Centre for Astrophysics, Alan Turing Building, Department of Physics and Astronomy, School of Natural Sciences,

The University of Manchester, Oxford Road, Manchester M13 9PL, UK

Department of Physics and Astronomy, University of British Columbia, Vancouver, BC V6T 1Z1, Canada

Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

Department of Physics, Korea Advanced Institute of Science and Technology (KAIST), 291 Daehak-ro, Yuseong-gu,

Daejeon 34141, Republic of Korea

Department of Physics, University of Toronto, 60 St. George Street, Toronto, ON M5S 1A7, Canada

David A. Dunlap Department of Astronomy, University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada

Brookhaven National Laboratory, Upton NY 11973-5000, USA

Kavli Institute for Particle Astrophysics and Cosmology and Physics Department, Stanford University, Stanford, CA 94305, USA

Department of Physics, University of Miami, 1320 Campo Sano Avenue, Coral Gables, FL 33146, USA

Department of Astronomy, University of Maryland, College Park, MD 20742, USA

Received 14 June 2024

Accepted 1 October 2024

ABSTRACT

The CO Mapping Array Project (COMAP) Pathfinder is performing line intensity mapping of CO emission to trace the distribution

of unresolved galaxies at redshift

∼

3. We present an improved version of the COMAP data processing pipeline and apply it to

the first two Seasons of observations. This analysis improves on the COMAP Early Science (ES) results in several key aspects. On

the observational side, all second season scans were made in constant-elevation mode, after noting that the previous Lissajous scans

were associated with increased systematic errors; those scans accounted for 50% of the total Season 1 data volume. In addition, all

new observations were restricted to an elevation range of 35–65 degrees to minimize sidelobe ground pickup. On the data processing

side, more e

ective data cleaning in both the time and map domain allowed us to eliminate all data-driven power spectrum-based

cuts. This increases the overall data retention and reduces the risk of signal subtraction bias. However, due to the increased sensitivity,

two new pointing-correlated systematic errors have emerged, and we introduced a new map-domain PCA filter to suppress these

errors. Subtracting only five out of 256 PCA modes, we find that the standard deviation of the cleaned maps decreases by 67%

on large angular scales, and after applying this filter, the maps appear consistent with instrumental noise. Combining all of these

improvements, we find that each hour of raw Season 2 observations yields on average 3.2 times more cleaned data compared to the ES

analysis. Combining this with the increase in raw observational hours, the e

ective amount of data available for high-level analysis is

a factor of eight higher than in the ES analysis. The resulting maps have reached an uncertainty of 25–50

μ

K per voxel, providing by

far the strongest constraints on cosmological CO line emission published to date.

Key words.

methods: data analysis – methods: observational – galaxies: high-redshift – di

use radiation –

radio lines: galaxies

1. Introduction

Line intensity mapping (LIM) is an emerging observational tech-

nique in which the integrated spectral line emission from many

Corresponding author;

j.g.s.lunde@astro.uio.no

unresolved galaxies is mapped in 3D as a tracer of the cosmolog-

ical large-scale structure (e.g., Kovetz et al. 2017, 2019). It rep-

resents a promising and complementary cosmological probe to,

say, galaxy surveys and cosmic microwave background (CMB)

observations. In particular, LIM o

ers the potential to survey

vast cosmological volumes at high redshift in a manner that is

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https:

creativecommons.org

licenses

4.0),

which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

A335, page 1 of 22

Lunde, J. G. S., et al.: A&A, 691, A335 (2024)

sensitive to emission from the entire galaxy population, not just

the brightest objects, as is the case for high-redshift galaxy sur-

veys (Bernal & Kovetz 2022). The most widely studied emis-

sion line for LIM purposes is the 21-cm line of neutral hydrogen

(Loeb & Zaldarriaga 2004; Bandura et al. 2014; Santos et al.

2017), which is the most abundant element in the universe,

but other high-frequency emission lines also appear promis-

ing (Korngut et al. 2018; Pullen et al. 2023; Akeson et al. 2019;

Crites et al. 2014; CCAT-Prime Collaboration 2022; Vieira et al.

2020; Karkare et al. 2022; Ade et al. 2020), in particular due to

their di

erent and complementary physical origin as well as

lower levels of astrophysical confusion, Galactic foregrounds,

and radio frequency interference.

The CO Mapping Array Project (COMAP) represents one

example of such an alternative approach and uses CO as the

tracer species (see, e.g., Lidz et al. 2011; Pullen et al. 2013;

Breysse et al. 2014). The COMAP Pathfinder instrument con-

sists of a 19-feed

focal plane array observing at 26–34 GHz

(Lamb et al. 2022) deployed on a 10.4 m Cassegrain telescope.

This frequency range corresponds to a redshift of

∼

4–

4 for the CO(1–0) line, a period during the Epoch of Galaxy

Assembly (Li et al. 2016). The Pathfinder instrument started

observing in 2019, and COMAP has previously published results

from the first year of data, named Season 1, in the COMAP

Early Science (ES) publications (Cleary et al. 2022; Lamb et al.

2022; Foss et al. 2022; Ihle et al. 2022; Chung et al. 2022;

Rennie et al. 2022; Breysse et al. 2022; Dunne et al. 2024).

These ES results provided the tightest constraints on the CO

power spectrum in the clustering regime published to date. Since

the release of the ES results, the COMAP Pathfinder instrument

has continued to observe while also implementing many impor-

tant lessons learned from Season 1, both in terms of observing

strategy and data processing methodology. Combining the obser-

vations from both Seasons 1 and 2 and improving the data anal-

ysis procedure, the new results improve upon the ES analysis

by almost an order of magnitude in terms of power spectrum

sensitivity.

This paper is the first of the COMAP Season 2 paper series,

and here we present the low-level data analysis pipeline and

map-level results derived from the full COMAP dataset as of

the end of Season 2 (November 2023). This work builds on

the corresponding Season 1 e

ort as summarized by Ihle et al.

(2022). The Season 2 power spectrum and null-test results are

presented by Stutzer et al. (2024), while Chung et al. (2024) dis-

cuss their cosmological implications in terms of structure forma-

tion constraints. In parallel with the Season 2 CO observations,

the COMAP Pathfinder continues to survey the Galactic plane,

with the latest results focusing on the Lambda Orionis region

(Harper et al. 2024).

The remainder of this paper is structured as follows: In

Sect. 2 we summarize the changes made to the observational

strategy in Season 2 and provide an overview of the current sta-

tus of data collection and accumulated volume. In Sect. 3 we

summarize our time-ordered data (TOD) pipeline with a focus

on the changes since the ES analysis. In Sect. 4 we study the sta-

tistical properties of the spectral maps produced by this pipeline

while paying particular attention to our new map-domain princi-

pal component analysis (PCA) filtering and the systematic errors

that this filter is designed to mitigate. In Sect. 5 we present

the current data selection methodology and discuss the resulting

improvements in terms of data retention in the time, map, and

power spectrum domains. In Sect. 6 we present updated trans-

We refer to a full detector chain as a “feed.”

fer function estimates and discuss their generality with respect

to non-linear filtering. Finally, we summarize and conclude in

Sect. 7.

2. Data collection and observing strategy

This section first briefly summarizes the COMAP instrument

and low-level data collection, which is extensively explored in

Lamb et al. (2022), before exploring the changes made to the

telescope and observing strategy between Seasons 1 and 2.

2.1. Instrument overview

The COMAP Pathfinder consists of a 19-feed 26–34 GHz spec-

trometer focal plane array fielded on a 10.4 m Cassegrain tele-

scope located at the Owens Valley Radio Observatory (OVRO).

At the central observing frequency of 30 GHz, the telescope has

a beam full width at half maximum (FWHM) of 4

′

. The 8 GHz-

wide RF signal is shifted to 2–10 GHz, and then split into two

4 GHz bands. The signals from each band are passed on to a sep-

arate “ROACH-2” field-programmable gate array spectrometer,

which further separates the 4 GHz-wide bands into 2 GHz-wide

sidebands. The spectrometer outputs 1024 frequency channels

for each of the four sidebands, for a spectral resolution of

∼

2 MHz.

The telescope is also equipped with a vane of microwave-

absorbing material, which is temporarily moved into the field of

view of the entire feed array before and after each hour-long

“observation” to provide absolute calibration of the observed

signal in temperature units. Each observation is stored in a single

HDF5 file containing both the spectrometer output and various

housekeeping data, and these files are referred to as “Level 1”

data. Each observation of one of the three target fields is divided

into ten to fifteen numbered “scans,” during which the telescope

oscillates in azimuth at constant elevation, repointing ahead of

the field to start a new scan each time the field has drifted through

the scanning pattern.

2.2. Status of observations

Table 1 shows the raw on-sky integration time per season.

COMAP Season 1 included 5200 on-sky observation hours col-

lected from May 2019 to August 2020, while the second sea-

son included 12 300 hours collected between November 2020

and November 2023. In these publications, we present results

based on a total on-sky integration time of 17 500 hours, a 3.4-

fold increase compared to the ES publications.

Several changes were made to the data collection and observ-

ing strategy before and during Season 2. Most of these changes

came as direct responses to important lessons learned during

the Season 1 data analysis and the aim was to increase the net

mapping speed, although one was necessary due to mechanical

telescope issues during Season 2. Overall, these changes were

highly successful, and Season 2 has a much higher data retention

than Season 1, which we discuss in Sect. 5. The most important

changes in the Season 2 observing strategy are the following:

1. Observations were restricted to an elevation range of 35

◦

–65

◦

in order to reduce the impact of ground pick-up via the tele-

scope’s sidelobe response. As discussed by Ihle et al. (2022),

the gradient of the ground pickup changes quickly at both

lower and higher elevations, and the corresponding observa-

tions were therefore discarded in the Season 1 analysis; in Sea-

son 2 we avoid these problematic elevations altogether.

A335, page 2 of 22

Lunde, J. G. S., et al.: A&A, 691, A335 (2024)

Table 1.

Overview of COMAP observation season definitions.

Name

Dates

Observing hours

Notes

Season 1

2019–08

2020

5200

Data published in Early Science. Contains 50% Lissajous, 50% CES

Season 2a

2020–04

2022

7900

100% CES from this point forward

Season 2b

2022–11

2023

4400

After azimuth drive slowdown and sampler frequency change

Notes.

Season 2 was split into two sub-seasons, respectively denoted 2a and 2b, as the telescope scanning speed was significantly reduced in May

2022 for mechanical reasons.

2. Similarly, Lissajous scans were abandoned in favor of solely

using constant elevation scans (CES), since Foss et al. (2022)

found elevation-dependent systematic errors associated with

the former.

3. The two frequency detector sub-bands, that previously

covered disjoint ranges of 26–30 GHz and 30–34 GHz

(Lamb et al. 2022), were widened slightly, such that they

now overlap; this mitigates data loss due to aliasing near the

band edges.

4. The acceleration of the azimuth drive was halved to increase

the longevity of the drive mechanism, which started to show

evidence of mechanical wear.

The latter two changes were only implemented in the second half

of Season 2, and they mark the beginning of what we refer to as

Season 2b. These changes are now discussed in greater detail.

2.3. Restricting the elevation range

Sidelobe pickup of the ground is a potentially worrisome sys-

tematic error for COMAP, especially since it is likely to be

pointing-correlated. Even though ground pickup is primarily

correlated with pointing in alt-azimuthal coordinates, the daily

repeating pointing pattern of COMAP means there will still be

a strong correlation in equatorial coordinates. Analysis of Sea-

son 1 observations (Foss et al. 2022), that ranged from

∼

◦

∼

◦

elevation, showed pointing-correlated systematic errors at

the highest and lowest elevations.

To study this e

ect in greater detail, we developed a set

of antenna beam pattern simulations using GRASP

for the

COMAP telescope (Lamb et al. 2022), and these showed the

presence of four sidelobes resulting from the four secondary

support legs (SSL), with each sidelobe o

set by

∼

◦

from

the pointing center. These simulations were convolved with the

horizon elevation profile at the telescope site, and the results

from these calculations are shown in Fig. 1. This figure clearly

shows that Fields 2 and 3 experience a sharp change in ground

pickup around 65

◦

–70

◦

elevation, as one SSL sidelobe transi-

tions between ground and sky. At very low elevations the ground

contribution also varies rapidly for all fields as two of the other

SSL sidelobes approach the horizon. While the low-level TOD

pipeline removes the absolute signal o

set per scan, gradients

in the sidelobe pickup over the duration of a scan still lead to

signal contamination. We have therefore restricted our observa-

tions to the elevation range of 35

◦

–65

◦

, where one SSL sidelobe

remains pointed at the ground, and the other three SSL side-

lobes are safely pointing at the sky, leaving us with a nearly

constant ground pickup. This change incurred little loss in obser-

vational e

ffi

ciency, as almost all allocated observational time

outside the new range could be reallocated to other fields within

the range.

https://www.ticra.com/software/grasp/

2.4. Abandoning Lissajous scans

The first season of observations contained an even distribu-

tion of Lissajous and constant elevation scans (CES), with the

aim of exploring the strengths and weaknesses of each. The

main strength of the Lissajous scanning technique is that it

provides excellent cross-linking by observing each pixel from

many angles, which is useful for suppressing correlated noise

with a destriper or maximum likelihood mapmaker. The main

drawback of this observing mode is that the telescope elevation

varies during a single scan, resulting in varying atmosphere and

ground pick-up contributions. In contrast, the telescope elevation

remains fixed during a CES, producing a simpler pick-up contri-

bution although with somewhat worse cross-linking properties.

When analyzing the Season 1 power spectra resulting from

each of the two observing modes, Ihle et al. (2022) find that the

Lissajous data both produce a highly significant power spec-

trum, especially on larger scales and fail key null tests. The

CES scans, on the other hand, produce a power spectrum con-

sistent with zero, and pass the same null tests. We therefore con-

clude that the significance in the Lissajous power spectrum is

due to residual systematic errors. Additionally, the main advan-

tage of Lissajous scanning, namely better cross-linking, prove

virtually irrelevant because of a particular feature of the COMAP

instrument: because all frequency channels in a single COMAP

sideband are processed through the same backend, the instru-

mental 1

gain fluctuations are extremely tightly correlated

across each sideband. As a result, the low-level TOD pipeline

is capable of simultaneously removing virtually all correlated

noise from both gain and atmosphere by common-mode subtrac-

tion (see Sect. 3.5). At our current sensitivity levels, we therefore

find no need to employ a complex mapmaking algorithm that

fully exploits cross-linking observations, but we can rather use

a computationally faster binned mapmaker (Foss et al. 2022).

After Season 1 we therefore concluded that there was no strong

motivation to continue with Lissajous scans, and in Season 2 we

employ solely CES.

2.5. Widening of frequency bands for aliasing mitigation

As discussed in detail by Lamb et al. (2022), the COMAP instru-

ment exhibits a small but non-negligible level of signal aliasing

near the edge of each sideband. In the Season 1 analysis, this

was accounted for simply by excluding the channels with alias-

ing power from other channels suppressed by less than 15 dB.

In total, 8% of the total COMAP bandwidth is lost due to this

ect, and this leads to gaps in the middle of the COMAP fre-

quency range. Both the origin of the problem and its ultimate

solution were known before the Season 1 observations started

(Lamb et al. 2022), but this took time to implement.

Band-pass filters applied after the first downconversion and

low-pass filters applied after the second downconversion allow

significant power above the Nyquist frequency into the sampler.

A335, page 3 of 22

Lunde, J. G. S., et al.: A&A, 691, A335 (2024)

100

150

200

250

300

350

Azimuth [degrees]

Elevation [degrees]

Field 1

Field 2

Field 3

Ground pickup [Kelvin]

Fig. 1.

Approximate sidelobe ground pickup as a function of az and el

pointing, simulated by convolving a beam (simulated using GRASP)

with the horizon profile (shown in gray) at the telescope site. The paths

of the three fields across the sky, in half-hour intervals, are shown, as

well as the Season 2 elevation limits at 35

◦

–65

◦

. These limits ensure

minimal ground pickup gradient across the field paths.

This is then aliased into the 0–2.0 GHz observing baseband,

requiring the contaminated channels to be excised. By increas-

ing the sampling frequency from 4.000 GHz to 4.250 GHz, the

Nyquist frequency is raised to 2.125 GHz, closer to the filter

edges. Not only is the amount of aliased power reduced, it is also

folded into frequencies above the nominal width of each 2 GHz

observing band. The existing samplers were able to accommo-

date the higher clock speed, but the field-programmable gate

array (FPGA) code had to be carefully tuned to reliably process

the data. This was finally implemented from the start of Sea-

son 2b, and the aliased power is now shifted outside the nom-

inal range of each band, such that the a

ected channels can be

discarded without any loss in frequency coverage. The number

of channels across the total frequency range is still 4096, so

the “native” Season 2b channels are 2.075 MHz wide, up from

1.953 MHz.

2.6. Azimuth drive slowdown

It became clear during Season 2 that the performance of the tele-

scope’s azimuth drives had degraded, owing to wear and tear

on the drive mechanisms caused by the telescope’s high accel-

erations. In order to protect the drives from damage, the analog

acceleration limit was reduced until the stress was judged by its

audible signature to be acceptable. Though not carefully quanti-

fied, this was about an order of magnitude change, and the min-

imum time for a scan is therefore about a factor of three less.

Additionally, the maximum velocity was reduced by a factor of

two in the drive software.

Figure 2 illustrates the old (Season 2a) and new (Season 2b)

pointing patterns, with the new pattern being slightly wider and

around four times slower. The new realized pointing pattern is

now also closer to sinusoidal since the drives are better able to

’keep up’ with the sinusoidal pattern of the commanded position,

due to the slower velocity.

With the new actually sinusoidal scanning pattern, the inte-

gration time is less uniform across each field in each observa-

tion, as the telescope spends more time pointing near the edges

of the field than it previously did. However, co-adding across the

observing season does smooth out the uneven integration time,

based on the receiver field of view (each of the 19 feeds observes

the sky at a position that is o

set from the others) and field rota-

tion (the telescope observes the fields at di

erent angles as they

move across the sky).

171.0

170.5

170.0

169.5

Right Ascension [degrees]

52.0

52.5

53.0

Declination [degrees]

S2a scan

S2b scan

Fig. 2.

Comparison of the faster pointing pattern from a Season 2a scan,

and the slower pointing pattern of a Season 2b scan. Both patterns show

a 5.5-minute constant elevation scan, as the field drifts across.

2.7. Data storage

With 19 feeds, 4096 native frequency channels, and a sampling

rate of about 50 samples

s, COMAP collects 56 GB

hour of raw

24-bit integer data, stored losslessly as 32-bit floats. The full

set of these raw data (combined with telescope housekeeping),

named “Level 1”, currently span about 800 TB. These data are

then filtered by our TOD pipeline into the so-called Level 2 data

(Foss et al. 2022), in which a key step is to co-add the native

frequency channels to 31.125 MHz width. These downsampled

channels form the basis of the higher-level map-making and

power spectrum algorithms. The total amount of Level 2 data is

about 50 TB. Both Level 1 and Level 2 files are now losslessly

compressed using the GZIP algorithm (Gailly & Adler 2023),

that achieves average compression factors of 2.4 and 1.4, respec-

tively, reducing the e

ective sizes of the two datasets to 330 TB

and 35 TB. The lower compression factor of the filtered data

is expected because the filtering leaves the data much closer to

white noise, and therefore with a much higher entropy.

The files are stored in the HDF5 format (Koranne 2011),

which allows seamless integration of compression. Both com-

pression and decompression happen automatically when writing

to and reading from each file. GZIP is also a relatively fast com-

pression algorithm, taking roughly one hour to compress each

hour of COMAP data on a typical single CPU core. Decompres-

sion takes a few minutes per hour of data, which is an insignifi-

cant proportion of the total pipeline runtime. HDF5 also allows

for arbitrary chunking of datasets before compression. Chunking

aids in optimizing performance since the Level 1 files consist of

entire observations (1 hour), but the current TOD pipeline imple-

mentation (see Sect. 3.9) reads only individual scans of 5 min-

utes each. We partition the data into chunks of 1-minute inter-

vals to minimize redundant decompression when accessing sin-

gle scans. Other compression algorithms have been tested, and

some, such as lzma

, achieve up to a 20% higher compression

factor on our data. They are, however, also much slower at both

compression and decompression, and they interface less easily

with HDF5.

3. The COMAP TOD pipeline

This section lays out the COMAP time-domain pipeline, named

l2gen

, focusing on the changes from the first generation

https://tukaani.org/xz/

A335, page 4 of 22

Lunde, J. G. S., et al.: A&A, 691, A335 (2024)

level2

data

l2gen

accept_mod

tod2comap

maps

comap2fpsx

mapPCA

power

spectra

level1

data

cleaned

maps

Fig. 3.

Flowchart of the COMAP pipeline, from raw Level1 data to final power spectra. Data products are shown as darker boxes and pipeline code

as lighter arrows.

l2gen

performs the time-domain filtering, turning raw data into cleaned Level 2 files.

accept_mod

performs scan-level data

selection on cleaned data.

tod2comap

is a simple binned mapmaker.

mapPCA

performs a map-level PCA filtering. Finally,

comap2fpsx

calculates

power spectra as described by Stutzer et al. (2024).

pipeline, which is described in detail by Foss et al. (2022). The

pipeline has been entirely re-implemented (see Sect. 3.9) for per-

formance and maintainability reasons but remains mathemati-

cally similar. Figure 3 shows a flowchart of the entire COMAP

pipeline, of which

l2gen

is the first step.

The purpose of the TOD pipeline is to convert raw detec-

tor readout (Level 1 files) into calibrated time-domain data

(Level 2 files) while performing two key operations: sub-

stantially reducing correlated noise and systematic errors, and

calibrating to brightness temperature units. COMAP uses a filter-

and-bin pipeline, meaning that we perform as much data clean-

ing as possible in the time domain, before binning the data into

maps with naïve noise-weighting. This leaves us with a biased

pseudo-power spectrum, that can be corrected by estimating the

pipeline transfer function (Sect. 6).

The following sections explain the main filters in the TOD

pipeline. The normalization (Sec. 3.3), 1

filter (Sec. 3.5), cal-

ibration and downsampling (Sect. 3.8) steps remain unchanged

from the ES pipeline, and are briefly summarized for complete-

ness. We denote the data at di

erent stages of the pipeline as

name

ν,

, where the

ν,

subscript indicates data with both frequency

and time dependence.

3.1. System temperature calculation

The first step of the pipeline is to calculate the system temper-

ature

sys

of each channel in the TOD. At the beginning and

end of each observation, a calibration vane of known tempera-

ture is moved into the field of view of all feeds. The measured

power from this “hot load”,

hot

, and the temperature of the vane,

hot

, are interpolated between calibrations to the center of each

scan. Power from a “cold load”,

cold

, is calculated as the aver-

age power of individual sky scans. The system temperature is

then calculated as (see Foss et al. 2022 for details)

sys

hot

−

CMB

hot

cold

−

(1)

under the approximation that the ground, sky, and telescope

share the same temperature.

3.2. Pre-pipeline masking

In ES,

l2gen

performed all frequency channel masking toward

the end of the pipeline. While some masks are data-driven

(specifically driven by the filtered data), others are not. We now

apply the latter category of masks prior to the filters, to improve

the filtering e

ectiveness. These are

– masking of channels that have consistently been found to be

correlated with systematic errors, and have been manually

flagged to always be masked;

– for data gathered before May 2022, masking of channels

with significant aliased power, as outlined in Sect. 2.5;

Frequency [GHz]

100

sys

[Kelvin]

System temperature

Running median

Running median + 5 Kelvin

Fig. 4.

Example of

sys

spike masking by running median. The system

temperature is shown in blue, the running median in orange, and the 5 K

cut above the running median in red. The

sys

values which are cut are

shown as a dashed instead of solid line.

– masking of channels with system temperature spikes, as out-

lined below.

The system noise temperature,

sys

, for each feed’s receiver

chain has a series of spikes at specific frequencies, believed to

result from an interaction between the polarizers and the corru-

gated feed horn (Lamb et al. 2022). The spikes are known to be

associated with certain systematic errors, and the a

ected fre-

quency channels are therefore masked out. The ES analysis used

a static system temperature threshold of 85 K for masking, but

the new version of

l2gen

applies a 400 channel-wide running

median kernel to the data and masks all frequencies with a noise

temperature of more than 5 K above the median. We repeat the

running median fit and threshold operation once on the masked

data, to reduce the impact of the spikes on the fit. The second

iteration uses a kernel width of 150 channels. The final running

median and threshold are illustrated in Fig. 4. As the system tem-

perature can vary quite a lot across the 4 GHz range, this method

fits the spikes themselves more tightly, while avoiding cutting

away regions of elevated but spike-free system temperature.

The spike frequencies vary from feed to feed, and we are

therefore not left with gaps in the redshift coverage of the final

3D maps. On average, we mask 6% of all frequency channels this

way. However, because the a

ected channels are, by definition,

more noisy, this only results in a loss of 3% of the sensitivity.

3.3. Normalization

The first filtering step in the TOD pipeline is to normalize the

Level 1 data by dividing by a low-pass filtered version of the

data and subtracting 1. The filter can be written as

norm

ν,

raw

ν,

〈

raw

〉

−

〈

raw

〉

−

{

·F{

raw

ν,

}}

(

knee

)

−

(2)

where

is the Fourier transform and

−

is the inverse Fourier

transform, both performed on the time-dimension of the data,

A335, page 5 of 22

Lunde, J. G. S., et al.: A&A, 691, A335 (2024)

Time [minutes]

Frequency [GHz]

Raw level1 data

Time [minutes]

Frequency [GHz]

After normalization

Time [minutes]

Frequency [GHz]

After 1/f filter

4 × 10

2 × 10

0 × 10

2 × 10

4 × 10

Power [digital units]

0.0010

0.0005

0.0000

0.0005

0.0010

Amplitude [normalized units]

0.0010

0.0005

0.0000

0.0005

0.0010

Amplitude [normalized units]

Fig. 5.

Illustration of the e

ect of TOD normalization and 1

filtering for a single feed and scan. The raw Level 1 data (left) are dominated

by frequency-dependent gain variations, that correspond to the instrumental bandpass. After normalization (middle), the signal is dominated by

common-mode gain fluctuations. Finally, after the 1

filter is applied (right), the common-mode 1

contribution has been suppressed, and the

data are dominated by white noise. The horizontal gray stripes indicate channels that were masked by the pipeline. All three stages happen before

absolute calibration, and the amplitudes are therefore given in arbitrary units.

and

is a low pass filter in the Fourier domain, with a spectral

slope

−

4, and a knee frequency

knee

01 Hz. This has the

ect of removing all modes on 100 s timescales and longer. As

the COMAP telescope crosses the entire field in 5–20 seconds,

the normalization has minimal impact on the sky signal in the

scanning direction, but heavily suppresses the signal perpendic-

ular to the scanning direction, as the fields take 5–7 minutes to

drift across.

The filter is performed per frequency channel, and the pri-

mary purpose of the normalization is to remove the channel-

to-channel variations in the gain, making channels more

comparable. After applying this filter, the white noise level in

each channel becomes the same, and common-mode 1

gain

fluctuations become flat across frequency. As a secondary con-

sideration, the normalization also removes the most slowly vary-

ing atmospheric and gain fluctuations, on timescales greater than

100 s.

The implications of the filter become more obvious by stat-

ing an explicit data model. The raw detector signal can be

modeled by the radiometer equation as a product of gain and

brightness temperature,

raw

ν,

, where both terms can in

theory vary freely. In practice, the LNA gain can be decomposed

into a mean time-independent part, and a multiplicative fluctua-

tion around this mean, such that we get

ν,

δg

), where

δg

is a small, frequency-independent term (often referred to as

gain noise). Because all frequencies are processed by the

same LNAs, these fluctuations become common-mode.

We can similarly decompose the brightness temperature into

a mean and fluctuation part. However, we can make no assump-

tion about the frequency-dependence of the fluctuation term, and

therefore simply write

ν,

. Here,

is now the sys-

tem temperature, as described in Eq. (1), and

ν,

are compara-

tively small fluctuations around this temperature. Putting this all

together we get the data model

raw

ν,

δg

)(

)

δg

ν,

)

(3)

≈

δg

ν,

)

where we have used the assumption that both fluctuations terms

are small to approximate

δg

ν,

≈

Under this data model, the low-pass filtered data is assumed

to take the form

〈

raw

〉

, such that the normalization filter

has the e

ect of transforming the data into normalized fluctua-

tions in gain and temperature around zero. Inserting Eq. (3) into

the filter as defined by Eq. (2), we get the normalized data

norm

ν,

δg

ν,

)

ν,

−

δg

ν,

(4)

Technically, we defined

and

to have no time-dependence,

but in the context of this filter, it makes sense to consider them

the temperature and gain terms that fluctuate more slowly than

∼

100 s timescales, as these are the timescales subtracted by this

filter.

The e

ect of the filter can be seen in the first two panels of

Fig. 5, which shows the TOD of a single scan in 2D before and

after the normalization. Before the normalization, the frequency-

dependent gain dominates, and the time variations are invisible.

After normalization, the data in each channel fluctuates around

zero.

3.4. Azimuth filter

Next, we fit and subtract a linear function in azimuth, to reduce

the impact of pointing-correlated systematic errors, first and

foremost being ground pickup by the telescope sidelobes. This

filter can be written as

ν,

norm

ν,

−

(5)

where

is fitted to the data per frequency, and az

is the azimuth

pointing of the telescope. Unlike in the ES pipeline, this filter

is now fitted independently for when the telescope is moving

eastward and westward, to mitigate some directional systematic

ects we have seen.

In Season 1, we also employed Lissajous scans, meaning that

an elevation term was also present in this equation. As we now

only observe in constant elevation mode, this term falls away.

3.5.

gain fluctuation filter

After normalization, the data are dominated primarily by gain,

and secondarily by atmospheric fluctuations, and both are

strongly correlated on longer timescales. Although the normal-

ization suppresses power on all timescales longer than 100 sec-

onds, we observe that common-mode noise still dominates the

total noise budget down to

∼

1 s timescales.

To suppress this correlated noise, we apply a specific 1

filter

by exploiting the simple frequency behavior of the gain

The filter is referred to as the polynomial filter in our ES publications.

A335, page 6 of 22

Lunde, J. G. S., et al.: A&A, 691, A335 (2024)

and atmosphere fluctuations. After we have normalized the data,

the amplitude of the gain fluctuations is the same across all fre-

quency channels, although fluctuating in time, as we see from the

term

δg

in Eq. (4). The atmospheric fluctuations will enter the

ν,

term of Eq. (4). To model the atmosphere, we approximate

both the atmosphere and the system temperature

as linear

functions of frequency, such that the combined term

ν,

itself linear. To jointly remove the gain and atmosphere fluctua-

tions, we therefore fitted and subtracted a first-order polynomial

across frequency for every time step:

ν,

point

ν,

−

(

)

(6)

where

and

are coe

ffi

cients fitted to the data each time step.

To ensure that both the atmosphere and the system temperature

are reasonably well approximated as linear in frequency, we fit

a separate linear polynomial to the 1024 channels of each of the

four 2 GHz sidebands.

This simple technique is remarkably e

ffi

cient at removing

noise, and we observe that the correlated noise is suppressed

by several orders of magnitude, after which white noise domi-

nates the uncertainty budget. This is illustrated in the last two

panels of Fig. 5. After the normalization, the signal is completely

dominated by common-mode 1

noise. After the 1

filter, the

correlated noise is e

ectively suppressed, and we are left with

almost pure white noise, as can be seen in the right panel, and

further discussed in Sect. 3.10.

3.6. PCA filtering

Principal component analysis (PCA) is a common and powerful

technique for dimensionality reduction (Pearson 1901). Given

a data matrix

ν,

, PCA produces an ordered basis

for the

columns of

ν,

, called the principal components of

ν,

. The

component amplitude can then be calculated by re-projecting the

components into the matrix, as

ν,

. For our purposes,

ν,

is the TOD, with frequencies as rows, and time-samples as

columns. The ordering of the principal components

is such

that the earlier components capture as much of the variance in

the columns of

ν,

as possible, and for any selected number of

components

comp

, the following expression is minimized:

∑

ν,

(

ν,

−

comp

∑

)

(7)

In other words, PCA provides a compressed version of

ν,

that approximates

ν,

as the sum of an ordered set of outer

products

ν,

≈

comp

∑

(8)

PCA is often employed on a dataset where the rows are inter-

preted as di

erent observations, and the columns are the multi-

dimensional features of these data. However, this is not a natural

interpretation for our purposes, and it makes more sense to sim-

ply look at PCA as a way of compressing a 2D matrix as a sum of

PCA has several equivalent interpretations and ways of solving for

the principal components. The principal components are, among other

things, the eigenvectors of the covariance matrix of

ν,

. This is how we

introduced the PCA in our ES publications. It is, however, both a slow

way of solving for the PCA components in practice and not the best

interpretation for our purposes.

Frequency [GHz]

0.00

0.05

Feed 6 freq. component

Time [minutes]

0.00

0.02

Amplitude

All feeds time component

Amplitude

Fig. 6.

Most significant component and amplitude of the all-feed PCA

filter applied on scan 3354205. The bottom plot shows the first PCA

component

, which is common to all feeds. The right plot shows the

corresponding amplitude

for feed 6. The outer product of these two

plots is shown in the central image and is the quantity subtracted from

the TOD by the filter. As the filter is applied on normalized data, the

amplitudes are all unitless, and the colorbar limits are

−

, half

the range of those in Fig. 5.

outer products – we have no special distinction between columns

and rows, and could equivalently have solved for the PCA of the

transpose of

ν,

, which would swap

and

A PCA is often performed because one is interested in

keeping the leading components, as these contain much of

the information in the data. However, we subtract the leading

components, because many systematic errors naturally decom-

pose well into an outer product of a frequency vector and a time

vector, while the CO signal does not (and is very weak in a single

scan).

In practice, we solve for the leading principal components

using a singular value decomposition algorithm (Halko et al.

2011), and then calculate the amplitudes as stated above. The

comp

leading components are then subtracted from the TODs,

leaving us with the filtered data

PCA

ν,

−

comp

∑

(9)

The Season 2 COMAP pipeline employs two time-domain PCA

filters, one of which was present in ES. In the following subsec-

tion, we introduce both filters and then explain how to decide the

number of leading components,

comp

, to subtract from the data.

Figure 6 shows an example of a strong component

picked

up by this filter in the bottom panel, and the corresponding com-

ponent amplitude

for feed 6 in the right panel. The center

image shows the resulting full frequency- and time-dependent

outer product between the two functions, which is the quantity

subtracted from the TOD. The time-dependent component is in

this example strongly temporally correlated, and could be resid-

ual atmospheric fluctuations, that typically have 1

-behavior

with a steep frequency slope.

The process of calculating the principal components and sub-

tracting them from the data constitutes a non-linear operation on

the data. This has the advantage of being much more versatile

against systematic errors that are di

ffi

cult to model using linear

filters, but the disadvantage is a more complicated impact on the

CO signal itself. This is further discussed in Sect. 6.4, where our

analysis shows that the PCA filter behaves linearly with respect

to any su

ffi

ciently weak signal, and that, at the scan-level, all

A335, page 7 of 22

Lunde, J. G. S., et al.: A&A, 691, A335 (2024)

Frequency [GHz]

0.05

0.00

0.05

Feed 14 freq. component

Time [minutes]

0.02

0.00

0.02

Amplitude

Feed 14 time component

Amplitude

Fig. 7.

Most significant component and amplitude of the per-feed PCA

filter applied on scan 3354205. The bottom plot shows the first PCA

component

for feed 14, while the right plot shows the corresponding

amplitude

for feed 14. The other product of the two quantities are

shown in the center plot, with colorbar limits of

−

plausible CO models (Chung et al. 2021) are su

ffi

ciently weak

by several orders of magnitude.

3.6.1. The all-feed PCA filter

The all-feed PCA filter, which was also present in the ES

pipeline, collapses the 19 feeds onto the frequency axis of the

matrix, producing a data matrix

ν,

of the 1

-filtered data with

a shape of (

feed

freq

TOD

)

(19

4096

∼

20 000) for a

scan with

feed

feeds,

freq

frequency channels, and

TOD

time

samples. The PCA algorithm outlined above is then performed

on this matrix. Combining the feed and frequency dimensions

means that a feature in the data will primarily only be picked up

by the filter if it is common (in the TOD) across all 19 feeds.

This is primarily the case for any atmospheric contributions, and

potentially standing waves that originate from the optics com-

mon to all feeds. It is, however, certainly not the case for the CO

signal, which will be virtually una

ected by this filter.

3.6.2. The per-feed PCA filter

The new per-feed PCA filter has been implemented to combat

systematic errors that vary from feed to feed. This filter employs

the PCA algorithm outlined above on each individual feed and is

performed on the output of the all-feed PCA filter. Additionally,

we found that downsampling the data matrix (using inverse vari-

ance noise weighing) by a factor of 16 in the frequency direction

before performing the PCA increased its ability to pick up struc-

tures in the data. The resulting data matrix

ν,

gets the shape

(

freq

TOD

)

(256

∼

20 000) for each feed. The down-

sampling is only used when solving for the time-domain com-

ponents

, and the full data matrix is used when calculating the

frequency amplitudes,

Targeting each feed individually makes us more susceptible

to CO signal loss, but the low signal-to-noise ratio (S

N), com-

bined with the fact that the CO signal cannot be naturally decom-

posed into an outer product, makes the impact on the CO signal

itself minimal. This filter appears to primarily remove compo-

nents consistent with standing waves from the individual optics

and electronics of each feed.

Figure 7 shows a typical strong component picked up by the

per-feed PCA filter for feed 14, similar to what was shown for

Sing. value number

1.0

1.1

1.2

Sing. value /

All-feed PCA filter

Per-feed PCA filter

Singular value cutoff

Number of PCA components removed

20%

40%

60%

80%

Percentage of scans

All-feed PCA filter

Per-feed PCA filter

Fig. 8.

Top

: largest singular values of the all-feed and per-feed PCA

filters, divided by

, for a random selection of ten scans. PCA com-

ponents with relative values above one are removed from the data and

are marked with crosses in this plot.

Bottom

: number of PCA compo-

nents subtracted across all scans. At least two components are always

subtracted by the all-feed PCA filter.

the all-feed PCA filter in Fig. 6. The per-feed PCA filter typically

picks up somewhat weaker features than the all-feed PCA filter,

as it is applied after the all-feed PCA filter. The time-dependent

component in Figure 7 is noticably noisier than the equivalent

component Fig. 6, also owing to the fact that it is fit with the

data from 1 feed, as opposed to 19 for the all-feed filter. The

physical origin of the feature shown is unknown.

3.6.3. Number of components

In the ES pipeline, the number of PCA components was fixed at

four for the all-feed filter, and the per-feed filter did not exist. We

now dynamically determine the required number of components

for each filter, per scan. This allows us to use more components

when needed, removing more systematic errors, and fewer when

not needed, incurring a smaller loss of CO signal.

We subtract principal components until the components are

indistinguishable from white noise, that can be inferred from

the singular values of each component. Let

be the expectation

value of the largest singular value of a (

) Gaussian noise

matrix (see Appendix A for how this value is derived). We sub-

tract principal components until we reach a singular value below

. However, for safe measure, we always subtract a minimum of

2 components for the all-feed PCA filter, as the signal impact of

this filter is minimal.

Figure 8 shows typical singular values, relative to

, for a

random selection of scans, and a histogram of the number of

components employed across all scans. The average number of

PCA components subtracted is 2.3 and 0.5 for the all-feed and

per-feed PCA, respectively; the most common number of com-

ponents subtracted is the minimum allowed in each case: two and

zero. The top part of the figure also demonstrates that there is a

sharp transition between the singular values of the components

that actually pick up meaningful features from the signal, and the

remaining noise components. For reference, the quite significant

components shown in Fig. 6 and Fig. 7 had singular value to

ratios of 2.9 and 1.6, respectively.

A335, page 8 of 22

Lunde, J. G. S., et al.: A&A, 691, A335 (2024)

3.7. Data-inferred frequency masking

After the PCA filters, we perform dynamic masking of frequency

channels identified from the filtered TOD. This mainly consists

of masking groups of channels that have substantially higher cor-

relations between each other than expected from white noise,

explained in more detail in Foss et al. (2022). This is assumed to

be caused by substantial residuals of gain fluctuations or atmo-

spheric signal. We did this by calculating the correlation matrix

between the frequency channels and perform

tests on the

white noise consistency of boxes and stripes of various sizes

across this matrix. We also mask individual channels with a

standard deviation significantly higher than expected from the

radiometer equation.

After the frequency masking, the 1

filter and PCA filters

are reapplied to the data, to ensure that their performance was

not degraded by misbehaving channels. The normalization and

pointing filters do not need to be reapplied, as they work inde-

pendently on each frequency channel.

3.8. Calibration and downsampling

The final step of the pipeline is to calibrate the data to tempera-

ture units and decrease the frequency resolution. After the nor-

malization step, the data are in arbitrary normalized units. Using

the system temperature calculated in Sect. 3.1 we calibrate each

channel of the data,

cal

ν,

sys

PCA

ν,

(10)

Finally, we downsample the frequency channels from 4096

native channels to 256 science channels. The Seasons 1 and

2a frequency channels are 1

953 MHz wide, while they are

075 MHz wide in Season 2b, after the change in sampling fre-

quency (Sect. 2.5). In both cases, the channels are downsampled

to a grid of 31

25 MHz, that exactly corresponds to a factor 16

downsampling for the older data. For the newer data, either 15

or 16 native channels will contribute to each science channel,

decided by their center frequency. The downsampling is per-

formed with inverse-variance weighting, using the system tem-

perature as uncertainty.

3.9. Implementation and performance

While the ES pipeline was written in Fortran, the Season 2

pipeline has been rewritten from scratch to run in Python.

Performance-critical sections are either written in C

and exe-

cuted using the Ctypes package or employ optimized Python

packages such as SciPy. Overall, the serial performance is simi-

lar to the ES pipeline, but the Season 2 pipeline employs a more

fine-grained and optimal MPI

OpenMP parallelization, making

it much faster on systems without a very large memory-to-core

ratio.

The pipeline is run on a small local cluster of 16 E7-8870v3

CPUs, with a total of 288 cores, in about a week of wall time,

totaling around 40 000 CPU hours for the full COMAP dataset.

The time-domain processing dominates this runtime, with a typi-

cal scan taking around 20–25 minutes to process on a single CPU

core.

3.10. Time-domain results

The Level 2 TODs outputted by

l2gen

are assumed to be almost

completely uncorrelated in both time and frequency dimensions,

such that the TOD are well approximated as white noise. To

PSD [normalized]

Raw real data (Level 1)

Filtered real data (Level 2)

White noise simulations

Slow scanning speed

Fast scanning speed

[Hz]

0.95

1.0

1.05

Fig. 9.

Average temporal power spectrum of unfiltered (green) and fil-

tered scans (blue), compared to correspondingly filtered white noise

simulations (red). The

-axis is broken at 1.06 and is logarithmic above

this. The data are averaged over scans, feeds, and frequencies, and are

normalized with respect to the highest

-bin. For context, the old and

new scanning frequencies (once across the field) are shown as vertical

lines (dot-dashed and dashed, respectively).

quantify the correlations in the time domain we calculate the

temporal power spectrum of each individual channel for all

scans. Figure 9 shows this power spectrum averaged over both

scans and frequencies, compared both to the equivalent power

spectrum of un-filtered Level 1 data, and to that of a TOD

obtained by injecting pure white noise in place of our real data

into the TOD pipeline.

Since the pipeline filtering removes more data on longer

timescales, the white noise simulation (red) gradually devi-

ates from a flat power spectrum on longer timescales, until it

falls rapidly at timescales below

∼

03 Hz due to the high-pass

normalization performed in the pipeline. The power spectrum

for the filtered real data (blue) also follows the same trend

on short timescales, but then increases on timescales around

2 Hz, due to small residual 1

gain fluctuations and

atmosphere remaining after the processing. Below

∼

03 Hz,

this power spectrum again falls rapidly due to the high-pass

filter. The power spectrum obtained from raw data that have

not gone through any filtering (green), simply increases on

longer timescales as expected, due to 1

gain and atmospheric

fluctuations.

The di

erence between the blue and red spectra shows the

residual correlated noise left in the data. While there is some

residual correlated noise, it is an insignificant fraction of our final

noise budget. We find that our real filtered data have a standard

deviation only 1.7% higher than that of the filtered white noise.

Compared to the amount of power we see in the unfiltered data,

we see that our pipeline is very e

ffi

cient at suppressing correlated

noise.

4. Mapmaking and map domain filtering

4.1. The COMAP mapmaker

COMAP employs a simple binned inverse-variance noise-

weighted mapmaker, identical to the one in ES (Foss et al. 2022).

This can be written as

ν,

∑

∈

ν,

/σ

ν,

∑

∈

/σ

ν,

(11)

A335, page 9 of 22

Lunde, J. G. S., et al.: A&A, 691, A335 (2024)

where

ν,

is an individual map voxel

ν,

represents the time-

domain data over time samples,

ν,

is the time-domain white

noise uncertainty, and

∈

means all time-samples

that

observes pixel

. We assume that the white noise uncertainty

ν,

is constant (for a single feed and frequency) over the duration of

a scan, and calculate it per scan as

√

Var(

ν,

−

ν,

−

)

(12)

then let

ν,

for all time-samples within the scan. This value

is also binned into maps, and is used as the uncertainty estimate

of the maps throughout the rest of the analysis:

ν,

∑

∈

/σ

(13)

In practice, we calculate per-feed maps, both because the map-

PCA filter (introduced in the next section) is performed on per-

feed maps, and because the cross-spectrum algorithm utilizes

groups of feed-maps.

The reason for not using more sophisticated mapmaking

schemes, such as destriping (Keihänen et al. 2010) or maximum

likelihood mapmaking (Tegmark 1997), is partially of necessity

– the COMAP TOD dataset is many hundreds of TB, making

iterative algorithms di

ffi

cult. However, the TOD pipeline has

also proven remarkably capable of cleaning most unwanted sys-

tematic errors from the data, especially correlated 1

noise,

as we saw from Fig. 9. As the main purpose of more sophis-

ticated mapmaking techniques is dealing with correlated noise,

COMAP is served well with a simple binned mapmaker.

The mapmaking algorithm is identical to the one used for

the ES analysis. However, as with

l2gen

the actual implementa-

tion has been rewritten from scratch in Python and C

, with a

focus on optimal parallelization and utilization of both MPI and

OpenMP.

4.2. Map-domain PCA filtering

The pipeline now employs a PCA filtering step also in the map

domain, in addition to the one we apply at the TOD level.

This technique is almost entirely analogous to the PCA fore-

ground subtraction often employed in 21 cm LIM experiments

(Chang et al. 2010; Masui et al. 2013; Anderson et al. 2018),

although we did not employ it to subtract foregrounds. The pri-

mary purpose of this filter is to mitigate a couple of pointing-

correlated systematic errors (see the next subsection) which

proved challenging to remove entirely in the time domain. The

method is similar to the TOD PCA algorithm from Sect. 3.6,

but instead of having the TOD data matrix

ν,

, we have a map

ν,

with one frequency and one (flattened) pixel dimension. The

data matrix then gets the shape (

)

(256

14 400), with

256 being the number of frequency channels in the map

and

120

14 400 the number of pixels in each

frequency slice (although many of the pixels in each individual

feed-map are never observed).

The technique we employ here is technically a slight gen-

eralization of the PCA problem, as we want to weigh individ-

ual voxels by their uncertainty when solving for the components

and amplitudes

. This is not possible in the regular PCA frame-

A voxel here is the 3D equivalent of a pixel, with two angular dimen-

sions and a redshift (frequency) dimension. Here we separate the voxel

dimensions into frequencies

and pixels

For the time-domain PCA, it was enough to weight individual chan-

nels, with all time-samples in that channel sharing the same weight.

This can be done with a normal PCA, as we show in Appendix B.

work without also morphing the modes one is trying to fit, as we

explain in Appendix B. As shown in Sect. 3.6, the first principal

component

and its amplitude

are the vectors that minimize

the value of the expression

∑

ν,

(

ν,

−

)

(14)

In other words, they are the two vectors such that their outer

product explains as much of the variance in

ν,

as possible.

This formulation of the PCA makes it obvious how to generalize

the problem to include weighting for individual matrix elements:

we can simply minimize the following sum,

∑

ν,

(

ν,

−

)

ν,

(15)

where

ν,

is the uncertainty in each voxel. Minimizing Eq. (15)

gives us the vectors

and

for which the outer product

explains as much of the variance in

ν,

as possible, weighted by

ν,

. The resulting map

ν,

−

represents the filtered map.

The process can then be repeated any number of times, solv-

ing for and subtracting a new set of vectors. We minimize the

expression in Eq. (15) iteratively with an algorithm outlined in

Appendix B, where we also explain why this is not equivalent to

simply performing the PCA on a noise-weighted map

ν,

/σ

ν,

Due to the large similarity of our technique to a regular PCA, we

simply refer to this filter as a PCA filter.

As for the TOD PCA filter, a selected number of components

are subtracted from the data maps

mPCA

ν,

−

comp

∑

(16)

This filtering is performed per feed, as the systematic errors

outlined in the next subsection manifest di

erently in di

erent

feeds. We have chosen

comp

5, which is further explained in

Sect. 4.3.3. Because the COMAP scanning strategy stayed the

same throughout Seasons 1 and 2a but changed with the azimuth

slowdown of Season 2b, we apply the map-PCA separately to

the former and latter, as the pointing-correlated e

ects we are

trying to remove might also be di

erent.

4.3. Newly discovered systematic effects

The two most prominent new systematic errors discovered in

the second season of observations have been dubbed the “turn-

around” and “start-of-scan” e

ects. They have in common that

they are di

ffi

cult to model in the time domain, subtle in individ-

ual scans, but strongly pointing-correlated, and they therefore

show up as large-scale features in the final maps. Additionally,

they are present to varying extents in all feeds, have similar quan-

titative behavior in the map-domain, and can both be removed

ectively with the map-PCA. The e

ects are discussed in the

subsections below, with further analysis shown in Appendix C.

4.3.1. The turn-around effect

The so-called “turn-around” e

ect can be observed as strongly

coherent excess power near the edges of the scan pattern, where

the telescope reverses direction in azimuth. Illustrations of this

ect can be seen in the first row, and partially in the third row,

of Fig. 10. The feature manifests at the top and bottom of the

A335, page 10 of 22

Lunde, J. G. S., et al.: A&A, 691, A335 (2024)

171

170

169

RA [degrees]

DEC [degrees]

before map-PCA

Feed 9, 26.89 GHz

171

170

169

RA [degrees]

after map-PCA

Significance [ ]

171

170

169

RA [degrees]

DEC [degrees]

before map-PCA

Feed 15, 31.27 GHz

171

170

169

RA [degrees]

after map-PCA

Significance [ ]

171

170

169

RA [degrees]

DEC [degrees]

before map-PCA

Feed 6, 27.11 GHz

171

170

169

RA [degrees]

after map-PCA

Significance [ ]

171

170

169

RA [degrees]

DEC [degrees]

before map-PCA

Feed 4, 33.58 GHz

171

170

169

RA [degrees]

after map-PCA

Significance [ ]

Fig. 10.

Selection of individual frequency maps from Field 2 before

(

left

) and after (

right

) the map-domain PCA filter. Because of the

uneven sensitivity among the pixels, and to emphasize the relevant sys-

tematic e

ects, all maps have been divided by their white noise uncer-

tainty. From top to bottom, each row shows maps that are dominated by

(1) the turn-around e

ect; (2) the start-of-scan e

ect; (3) both e

ects

simultaneously; and (4) neither e

ect. Both e

ects appear to manifest

twice, on two slightly o

set maps. This o

set e

ect originates from the

physical placement of the feeds, as they observe the fields as they are

both rising and setting on the sky. In the equatorial coordinate system of

these maps, the telescope scans vertically, and the field drifts from left

to right.

maps, as this is where the telescope turn-around happens for

Field 2 in equatorial coordinates. The feature oscillates slowly

across the frequency domain, and the leading theory of its origin

is some standing wave oscillation induced by mechanical vibra-

tions. The e

ect is somewhat less pronounced in Season 2b, after

the reduction in telescope pointing speed and acceleration, but it

is still present.

4.3.2. The start-of-scan effect

A related e

ect is called the “start-of-scan” e

ect, that is a wave-

like feature in frequency that occurs at the beginning of every

scan and decays exponentially with a mean lifetime of around

19 s. As the telescope always starts each scan at the same side of

each sky field, this systematic e

ect shows up in the map domain

PCA component 1

PCA amplitude 1

PCA component 2

PCA amplitude 2

PCA component 3

PCA amplitude 3

PCA component 4

PCA amplitude 4

PCA component 5

Frequency [GHz]

PCA amplitude 5

Fig. 11.

Leading PCA components

(

left

) and their respective fre-

quency amplitudes

(

right

) for Field 2 as observed by feed 6 prior to

the map-PCA filter; this feed is the most sensitive to pointing-correlated

systematics. All maps are divided by their respective uncertainties to

highlight the key morphology. All rows share the same color range and

-axis scale, but the specific values have been omitted as they are not

easily interpretable.

as a strong feature on the Eastern edge of the map, as can be

seen in the second and third rows of Fig. 10. Next to the strong

positive or negative signal (this varies by frequency) at the very

edge of the map, the opposite power, at lower amplitude, can be

observed as we move Westward across the map. This opposite

power is simply a ringing feature from the normalization per-

formed during the pipeline (see Appendix C for details).

The exact origins of the “start-of-scan” e

ect are unknown,

but the fact that it only happens at the beginning of scans (that are

separated by a repointing to catch up with the field), and disap-

pears during constant elevation scanning, suggests that a poten-

tial candidate is mechanical vibrations induced by the repointing.

We also observe the e

ect to be mostly associated with one of the

four DCM1s (first downconversion module), namely DCM1-2,

relating to feeds 6, 14, 15, 16, and 17. The e

ect’s strong corre-

lation with DCM1-2 points to a possible source in the local oscil-

lator cable, shared by the channels in a DCM1 module; imper-

fect isolation of the mixer would cause a weak common-mode

resonance to manifest.

An important detail in this analysis is that, because of the

normalization and 1

filter in the TOD pipeline, any stand-

ing wave signal with a constant resonant cavity wavelength over

time will be filtered away. For a standing wave to survive the fil-

tering, it must have a changing wavelength. Prime suspects for

the origin of this e

ect are therefore optical cavities that could

expand or contract in size, or cables that could be stretched.

A335, page 11 of 22