# LiteBIRD Science Goals and Forecasts: constraining isotropic cosmic birefringence

E de la Hoz, P Diego-Palazuelos, J Errard, A Gruppuso, B Jost, R.M

Sullivan, M Bortolami, Y Chinone, L.T Hergt, E Komatsu, et al.

► To cite this version:

E de la Hoz, P Diego-Palazuelos, J Errard, A Gruppuso, B Jost, et al. LiteBIRD Science Goals and Forecasts: constraining isotropic cosmic birefringence. Journal of Cosmology and Astroparticle Physics, 2025, 07, pp.083. 10.1088/1475-7516/2025/07/083. hal-05029509

HAL Id: hal-05029509

https://hal.science/hal-05029509v1

Submitted on 16 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

LiteBIRD science goals and forecasts:
constraining isotropic cosmic birefringence

To cite this article: E. de la Hoz et al JCAP07(2025)083

View the [article online](#) for updates and enhancements.

You may also like

- Requirements on the gain calibration for LiteBIRD polarisation data with blind component separation
F. Carralot, A. Carones, N. Krachmalnicoff et al.

- Impact of beam far side-lobe knowledge in the presence of foregrounds for LiteBIRD
C. Leloup, G. Patanchon, J. Errard et al.

- Polarization angle requirements for CMB B-mode experiments. Application to the LiteBIRD satellite
The LiteBIRD collaboration, P. Vielva, E. Martínez-González et al.
# LiteBIRD science goals and forecasts: constraining isotropic cosmic birefringence

<figure><img src="image_2.png" /></figure>

## The LiteBIRD collaboration

E. de la Hoz,ª,bk P. Diego-Palazuelos,ᵇ J. Errard,ᶜ A. Gruppuso,ᵈ,ᵉ B. Jost,ᵉ,f,bl

R.M. Sullivan,ᵈ,g,h M. Bortolami,*,j Y. Chinone,*,f L.T. Hergt,*,l E. Komatsu,*,b,f

Y. Minami,⁹ I. Obata,f D. Paoletti,*,d,e D. Scott,г P. Vielva,⁹ D. Adak,⁰ R. Akizawa,⁹

A. Anand,q J. Aumont,*,r C. Baccigalupi,*,s,t,u A.J. Banday,*,r R.B. Barreiro,⁹ N. Bartolo,*,v,w,x

S. Basak,г A. Basyrov,г M. Bersanelli,*,z,aa T. Brinckmann,*,f F. Cacciotti,ab E. Calabrese,ac

P. Campeti,*,j,b,ad E. Carinos,*,r A. Carones,*,s,t F. Carralot,*,s,ae,t F. J. Casas,⁹ M. Citran,ᶜ

L. Clermont,af F. Columbro,ab,ag G. Coppi,ah A. Coppolecchia,ab,ag F. Cuttaia,⁴

P. de Bernardis,ab,ag M. De Lucia,ai M. De Petris,ab,ag S. Della Torre,aj E. Di Giorgi,ai

H.K. Eriksen,г E. Ferreira,г F. Finelli,*,d,e C. Franceschet,*,z,aa U. Fuskeland,г G. Galloni,*,i,q

M. Galloway,г M. Gerbino,*,j M. Gervasi,*,ah,aj R.T. Génova-Santos,*,o,ak T. Ghigna,*,k,f

S. Giardiello,ac C. Gimeno-Amo,ⁿ E. Gjerløw,г M. Hazumi,al,am,f,an S. Henrot-Versillé,¹

E. Hivon,ao H. Ishino,ap K. Kohri,al L. Lamagna,ab,ag M. Lattanzi,*,j C. Leloup,*,f

M. Lembo,*,f F. Levrier,*,q M. López-Caniego,*,ar,as G. Luzzi,*,at E. Martínez-González,ⁿ

S. Masi,ab,ag S. Matarrese,*,v,w,x,au S. Micheli,ab M. Migliaccio,*,q,av M. Monelli,*,f

L. Montier,*,r G. Morgante,*,d R. Nagata,am T. Namikawa,*,f P. Natoli,*,i,j A. Occhiuzzi,ab

L. Pagano,*,i,j,aw A. Paiella,ab,ag G. Pascual-Cisneros,*,f,n V. Pavlidou,*,ax,ay V. Pelgrims,az

F. Piacentini,ab,ag G. Piccirilli,г G. Polenta,*,at L. Porcelli,ba N. Raffuzzi,*,i M. Remazeilles,ⁿ

A. Ritacco,bb A. Rizzieri,*,bc,c J.A. Rubiño-Martín,*,o,ak M. Ruiz-Granda,*,n,bd Y. Sakurai,*,be,f

J. Sanghavi,*,bf,bg M. Shiraishi,be S.L. Stever,*,aw,ap,*,f Y. Takase,ap K. Tassis,*,ax,ay

L. Terenzi,*,d M. Tomasi,*,z,aa M. Tristram,*,l L. Vacher,*,s B. van Tent,*,l D. Watts,*,h

I.K. Wehus,*,h G. Weymann-Despres,*,bc,l B. Winter,*,bh,bi E.J. Wollack,*,bj and Y. Zhou,*,k

<figure><img src="image_3.png" /></figure>

<footnote>ªCentre Pierre Binetruy, CNRS-University of California, Berkeley, IRL2007 CPB-IN2P3, Berkeley, CA, U.S.A.</footnote>

<footnote>ᵇMax Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany</footnote>

<footnote>ᶜUniversité Paris Cité, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France</footnote>

<footnote>ᵈINAF — OAS Bologna, via Piero Gobetti, 93/3, 40129 Bologna, Italy</footnote>

<footnote>ᵉINFN Sezione di Bologna, Viale C. Berti Pichat, 6/2, 40127 Bologna, Italy</footnote>

<footnote>fKavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan</footnote>
q Department of Physics and Astronomy, University of British Columbia,
6224 Agricultural Road, Vancouver, BC V6T1Z1, Canada

h Institute of Theoretical Astrophysics, University of Oslo, Blindern, Oslo, Norway

i Dipartimento di Fisica e Scienze della Terra, Università di Ferrara,
Via Saragat 1, 44122 Ferrara, Italy

j INFN Sezione di Ferrara, Via Saragat 1, 44122 Ferrara, Italy

k International Center for Quantum-field Measurement Systems for Studies of the Universe and Particles (QUP),
High Energy Accelerator Research Organization (KEK),
Tsukuba, Ibaraki 305-0801, Japan

l Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France

m Research Center for Nuclear Physics, Osaka University, Ibaraki, Osaka, 567-0047, Japan

n Instituto de Fisica de Cantabria (IFCA, CSIC-UC),
Avenida los Castros SN, 39005, Santander, Spain

o Instituto de Astrofísica de Canarias, E-38200 La Laguna, Tenerife, Canary Islands, Spain

p The University of Tokyo, Department of Physics, Tokyo 113-0033, Japan

q Dipartimento di Fisica, Università di Roma Tor Vergata,
Via della Ricerca Scientifica, 1, 00133, Roma, Italy

r IRAP, Université de Toulouse, CNRS, CNES, UPS, Toulouse, France

s International School for Advanced Studies (SISSA), Via Bonomea 265, 34136, Trieste, Italy

t INFN Sezione di Trieste, via Valerio 2, 34127 Trieste, Italy

u IFPU, Via Beirut, 2, 34151 Grignano, Trieste, Italy

v Dipartimento di Fisica e Astronomia “G. Galilei”, Università degli Studi di Padova,
via Marzolo 8, I-35131 Padova, Italy

w INFN Sezione di Padova, via Marzolo 8, I-35131, Padova, Italy

x INAF, Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122, Padova, Italy

y School of Physics, Indian Institute of Science Education and Research Thiruvananthapuram,
Maruthamala PO, Vithura, Thiruvananthapuram 695551, Kerala, India

z Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, 20133, Milano, Italy

aa INFN Sezione di Milano, Via Celoria 16, 20133, Milano, Italy

ab Dipartimento di Fisica, Università La Sapienza, P. le A. Moro 2, Roma, Italy

ac School of Physics and Astronomy, Cardiff University, Cardiff CF24 3AA, U.K.

ad Excellence Cluster ORIGINS, Boltzmannstr. 2, 85748 Garching, Germany

ae Università di Trento, Dipartimento di Fisica, Via Sommarive 14, 38123, Trento, Italy

af Centre Spatial de Liège, Université de Liège, Avenue du Pré-Aily, 4031 Angleur, Belgium

ag INFN Sezione di Roma, P.le A. Moro 2, 00185 Roma, Italy

ah University of Milano Bicocca, Physics Department, p.zza della Scienza, 3, 20126 Milan, Italy

ai INFN Sezione di Pisa, Largo Bruno Pontecorvo 3, 56127 Pisa, Italy

aj INFN Sezione Milano Bicocca, Piazza della Scienza, 3, 20126 Milano, Italy

ak Departamento de Astrofísica, Universidad de La Laguna (ULL),
E-38206, La Laguna, Tenerife, Spain

al Institute of Particle and Nuclear Studies (IPNS),
High Energy Accelerator Research Organization (KEK),
Tsukuba, Ibaraki 305-0801, Japan

am Japan Aerospace Exploration Agency (JAXA), Institute of Space and Astronautical Science (ISAS),
Sagamihara, Kanagawa 252-5210, Japan
an The Graduate University for Advanced Studies (SOKENDAI),
Miura District, Kanagawa 240-0115, Hayama, Japan

ao Institut d'Astrophysique de Paris, CNRS/Sorbonne Université, Paris, France

ap Okayama University, Department of Physics, Okayama 700-8530, Japan

aq Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS,
Sorbonne Université, Université de Paris,
75005 Paris, France

ar Aurora Technology for the European Space Agency,
Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada,
Madrid, Spain

as Universidad Europea de Madrid, 28670, Madrid, Spain

at Space Science Data Center, Italian Space Agency, via del Politecnico, 00133, Roma, Italy

au Gran Sasso Science Institute (GSSI), Viale F. Crispi 7, I-67100, L'Aquila, Italy

av INFN Sezione di Roma2, Università di Roma Tor Vergata,
via della Ricerca Scientifica, 1, 00133 Roma, Italy

aw Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405, Orsay, France

ax Institute of Astrophysics, Foundation for Research and Technology — Hellas,
Vasilika Vouton, GR-70013 Heraklion, Greece

ay Department of Physics and ITCP, University of Crete, GR-70013, Heraklion, Greece

az Université Libre de Bruxelles, Av. Franklin Roosevelt 50, 1050 Bruxelles, Belgium

ba Istituto Nazionale di Fisica Nucleare-Laboratori Nazionali di Frascati (INFN-LNF),
Via E. Fermi 40, 00044, Frascati, Italy

bb Université Grenoble Alpes, CNRS, LPSC-IN2P3, 53, avenue des Martyrs, 38000 Grenoble, France

bc Department of Physics, University of Oxford,
Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, U.K.

bd Departamento de Física Moderna, Universidad de Cantabria,
Avda. los Castros s/n, E-39005 Santander, Spain

be Suwa University of Science, Chino, Nagano 391-0292, Japan

bf Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians Universität München,
Scheinerstr. 1, 81679 München, Germany

bg GRAPPA, Institute for Theoretical Physics Amsterdam, University of Amsterdam,
Science Park 904, 1098 XH Amsterdam, The Netherlands

bh Mullard Space Science Laboratory, Department of Space and Climate Physics,
University College London,
Gower Street, WC1E 6BT, London, U.K.

bi Physics and Astronomy Department, University College London (UCL), London, U.K.

bj NASA Goddard Space Flight Center, Greenbelt, MD 20771, U.S.A.

bk Department of Physics, University of California, Berkeley,
366 Physics North, Berkeley, CA, U.S.A.

bl Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo,
Kashiwa, Chiba 277-8583, Japan

E-mail: alessandro.gruppuso@inaf.it, josquin@apc.in2p3.fr

**ABSTRACT:** Cosmic birefringence (CB) is the rotation of the photons' linear polarisation plane during propagation. Such an effect is a tracer of parity-violating extensions of standard electromagnetism and would probe the existence of a new cosmological field acting as dark matter
or dark energy. It has become customary to employ cosmic microwave background (CMB) polarised data to probe such a phenomenon. Recent analyses on *Planck* and *WMAP* data provide a hint of detection of the isotropic CB angle with an amplitude of around 0.3° at the level of 2.4 to 3.6σ. In this work, we explore the *LiteBIRD* capabilities in constraining such an effect, accounting for the impact of the more relevant systematic effects, namely foreground emission and instrumental polarisation angles. We build five semi-independent pipelines and test these against four different simulation sets with increasing complexity in terms of non-idealities.

All the pipelines are shown to be robust and capable of returning the expected values of the CB angle within statistical fluctuations for all the cases considered. We find that the uncertainties in the CB estimates increase with more complex simulations. However, the trend is less pronounced for pipelines that account for the instrumental polarisation angles. For the most complex case analysed, we find that *LiteBIRD* will be able to detect a CB angle of 0.3° with a statistical significance ranging from 5 to 13 σ, depending on the pipeline employed, where the latter uncertainty corresponds to a total error budget of the order of 0.02°.

KEYWORDS: CMBR experiments, CMBR polarisation, cosmological parameters from CMBR, CMBR theory

arXiv ePRINT: 2503.22322
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Simulations</td><td>3</td></tr><tr><td>3 Techniques employed</td><td>6</td></tr><tr><td>3.1 Pipelines without control of instrumental polarisation angles</td><td>6</td></tr><tr><td>3.2 Pipelines with control of instrumental polarisation angles</td><td>10</td></tr><tr><td>4 Analyses</td><td>18</td></tr><tr><td>4.1 Pipelines without control of instrumental polarisation angles</td><td>18</td></tr><tr><td>4.2 Pipelines with control of instrumental polarisation angles</td><td>23</td></tr><tr><td>5 Discussion about the impact of foreground contamination</td><td>38</td></tr><tr><td>6 Conclusion</td><td>40</td></tr><tr><td>A Definition of the template-based Minami-Komatsu estimator</td><td>42</td></tr></table>

## 1 Introduction

Parity-violating extensions of the standard electromagnetism due to the coupling between photons and a (pseudo-)scalar field, for instance, an axion, via a Chern-Simons term [1], can be traced with the cosmic birefringence (henceforth CB) effect, i.e. the in-vacuo<fnref n="1" /> rotation of the linear polarisation plane of photons during their propagation [2–6]. Such a rotation is expected to be null in the standard Maxwellian electromagnetism, so detecting a CB angle different from zero would probe the existence of a new cosmological scalar field, possibly acting as dark matter or dark energy [7, 8].

In particular, for ultra-light scalar fields, the Chern-Simons term predicts a CB rotation proportional to the slowly-rolling distance travelled by the photons [3, 9, 10]. Since cosmic microwave background (CMB) photons are linearly polarised at the last scattering surface due to Thomson scattering, and have travelled the furthest possible distance through the Universe, it has become customary to exploit polarised CMB data to probe for such an effect.

After COBE [11], WMAP [12], and Planck [13], LiteBIRD [14] is the next CMB space-based mission aiming to measure the polarised signal at large angular scales with unprecedented precision. This work is part of a series of papers presenting the science achievable by LiteBIRD space mission and expanding on the overview published in ref. [14]. In particular, we focus on CB and forecast LiteBIRD’s capabilities in constraining this phenomenon, taking into account realistic systematic effects of both instrumental and astrophysical origins.

Our present analysis is limited to the isotropic CB effect, where the rotation is constant over the whole sky and over time. Hence, we do not consider here either the isotropic CB effect with a time-varying rotation angle due to an oscillatory motion of the axion

<footnote n="1">¹Here, with “in-vacuo” we mean from the point of view of the standard model.</footnote>
background<fnref n="2" /> [19, 20], or the anisotropic CB effect, where the rotation angle depends on the direction of observation [21], and which has been found to be compatible with zero in previous studies [22–28]. In this case, the CMB angular power spectra are rotated by the isotropic CB angle $\beta$ in the following way, see e.g. refs. [29–31],

$$C_{\ell}^{TT,\text{obs}} = C_{\ell}^{TT}, \qquad (1.1)$$

$$C_{\ell}^{TE,\text{obs}} = C_{\ell}^{TE} \cos(2\beta), \qquad (1.2)$$

$$C_{\ell}^{TB,\text{obs}} = C_{\ell}^{TE} \sin(2\beta), \qquad (1.3)$$

$$C_{\ell}^{EE,\text{obs}} = C_{\ell}^{EE} \cos^2(2\beta) + C_{\ell}^{BB} \sin^2(2\beta), \qquad (1.4)$$

$$C_{\ell}^{BB,\text{obs}} = C_{\ell}^{BB} \cos^2(2\beta) + C_{\ell}^{EE} \sin^2(2\beta), \qquad (1.5)$$

$$C_{\ell}^{EB,\text{obs}} = \frac{1}{2} (C_{\ell}^{EE} - C_{\ell}^{BB}) \sin(4\beta), \qquad (1.6)$$

where the spectra with the “obs” label are observed and the spectra without the “obs” tag are those that would be observed in the case of null CB. Eqs. (1.3) and (1.6) ignore any pre-recombination TB and EB signals, as we consider CB the only parity-violating extension to the standard model (see e.g. ref. [32] for a review).

Unfortunately, CB investigations based solely on CMB data are limited due to a degeneracy between the instrument polarisation angles ($\alpha$) and the birefringence angle ($\beta$). This degeneracy is the main obstacle in performing such analyses, since the uncertainty in $\alpha$ typically dominates the total error budget for most of the currently available CMB observations [33, 34]. For example, the estimate of $\beta$ provided by the Planck collaboration [35] is compatible with the mentioned value of 0.3° but the instrumental uncertainty is of the same order of magnitude, whereas the statistical part is at least 6 times smaller. The degree of degeneracy depends on the mass of the axion [36–39]. However, as typically assumed in the literature and supported by current data [40], in our analysis we will consider such degeneracy to be complete. This would correspond to an axion mass small enough ($m_{\phi} \lesssim 10^{-32}$ eV) to make the axion oscillate only after the reionisation epoch, see e.g. refs. [39, 41].<fnref n="3" /> See also refs. [44, 45] for further theoretical interpretations of the EB Planck spectrum.

The most recent isotropic CB constraints, based on a novel technique [46–48] (hereafter the Minami-Komatsu or MK technique), also exploit the information contained in the foreground emission and simultaneously recover both the $\beta$ and $\alpha$ angles. This new method was proposed to break the degeneracy between $\beta$ and $\alpha$, finding hints of a detection of the isotropic CB angle $\beta \simeq 0.3^\circ$ [49–52] at a confidence level of 2.4 to 3.6$\sigma$. The MK technique represents a step forward in the investigation of parity-violating extensions of electromagnetism. Unfortunately, a deeper knowledge of the foregrounds themselves is necessary to provide a robust conclusion, and current results may yet be plagued by unknown systematics [50, 53–55]. Nevertheless, in ref. [56], using only reprocessed Planck-LFI and WMAP data, the authors still found $\beta$ compatible with the previously reported value of approximately 0.3°, primarily derived from

<footnote n="2">²Such dark-matter axions with mass scales around $10^{-21}$ eV $\lesssim m_{\phi} \lesssim 10^{-19}$ eV have been constrained by the polarisation oscillation measurements of BICEP/Keck [15, 16], SPT [17], and POLARBEAR [18].</footnote>

<footnote n="3">³See refs. [42, 43] for the dependence of the anisotropic birefringence effect on the axion mass also in the cross-spectra or cross-bispectra with CMB fields.</footnote>
Planck-HFI maps. It is interesting that this compatibility was found despite the very different astrophysical and instrumental systematics present in these data sets. Recently [57] found that their global polarization angle shows a 2.5 σ departure from zero. See also ref. [58] for a recent review of the latest constraints on the isotropic and anisotropic CB effects.

The paper is organised as follows. In section 2, we describe the simulations used, which contain realistic systematic effects coming both from astrophysical signals and from the instrument. We employ five different pipelines to estimate β, described in section 3.

In particular, we divide them into two categories: in section 3.1, we present the techniques that do not explicitly consider instrumental α, and, in section 3.2, we describe those that simultaneously estimate β and α. Section 4 presents the results from employing the different pipelines on the simulations and in section 5 we discuss the impact of foregrounds on our results. In section 6 we compare and discuss the results obtained, and finish with the main conclusions from this work.

## 2 Simulations

We consider four different sets of simulations, each with increasing complexity in terms of non-idealities.

We call these sets “Phases” and label them from 1 to 4. For each phase, we have produced 2200 polarised maps, i.e., 100 for each of the 22 LiteBIRD frequency channels. These are all provided at HEALPix<fnref n="4" /> [59] resolution $N_{side} = 512$ (which means that the number of pixels needed to cover a CMB map is $N_{pix} = 12N_{side}^2$). We convolved each frequency map with the Gaussian of its nominal beam width. See table 3 of ref. [14] for specifics of all the LiteBIRD frequency bands. Moreover, we assume that each frequency map is obtained with a Dirac delta function bandpass.

The simulated frequency maps $d_f$, in CMB units, are expressed as

$$ \mathbf{d}_f = \mathbf{R}^{\text{inst}}(\alpha_f)\mathbf{B}(\theta_{\text{beam}}^f) [\mathbf{R}^{\text{biref}}(\beta)\mathbf{A}_{\text{CMB}}, \mathbf{A}_{\text{dust}}^f, \mathbf{A}_{\text{sync}}^f] \begin{bmatrix} \mathbf{s}_{\text{CMB}} \\ \mathbf{s}_{\text{dust}} \\ \mathbf{s}_{\text{sync}} \end{bmatrix} + \mathbf{n}_f, \quad (2.1) $$

where boldface quantities are spanning the I, Q, and U Stokes parameters, and f labels the frequency band (in the rest of the text, we use the quantity ν for the frequency); $\mathbf{R}^{\text{inst}}(\alpha_f)$ is a rotation matrix modeling the imperfect $\alpha_f$ calibration, and $\mathbf{R}^{\text{biref}}(\beta)$ models the effect of CB. Note that $\mathbf{R}^{\text{biref}}$ is only applied to the CMB component $s_{\text{CMB}}$. For a given sky component, $\mathbf{R}^x$ is defined as:

$$ \mathbf{R}^{\mathbf{x}}(\gamma) = \begin{pmatrix} 1 & 0 & 0 \\ 0 \cos(2\gamma) & -\sin(2\gamma) & 0 \\ 0 \sin(2\gamma) & \cos(2\gamma) & 0 \end{pmatrix} \quad \text{with } \mathbf{x} \in \{\text{inst, biref}\}. \quad (2.2) $$

$\mathbf{B}$ is a 3 × 3 diagonal matrix beam operator that takes as input the beam size $\theta_{\text{beam}}^f$ and is applied to each frequency band and each Stokes parameter independently; $\mathbf{A}_{\text{CMB}}$, $\mathbf{A}_{\text{dust}}$ and $\mathbf{A}_{\text{sync}}$ contain the dust and synchrotron spectral energy densities (SEDs), parameterised

<footnote n="4">⁴http://healpix.sourceforge.net</footnote>
<table><caption>**Table 1.** Properties of the various simulations phases, in terms of miscalibration ($\alpha_f$) and birefringence ($\beta$) angles, and foreground spectral parameters ($\beta_d$, $T_d$, and $\beta_s$). In the case of Phase 3 and 4, the random $\alpha_f$ are drawn from uniform distributions $\alpha_f \in [-\sigma_{\alpha_f}, \sigma_{\alpha_f}]$ which widths are specified in table 2 and derived in ref. [60].</caption><thead><tr><th>Phase</th><th>α<sub>f</sub> [deg]</th><th>β [deg]</th><th>β<sub>d</sub>, T<sub>d</sub>, β<sub>s</sub></th></tr></thead><tbody><tr><td>1</td><td>0</td><td>0.0</td><td>Spatially constant (pysm d0s0)</td></tr><tr><td>2</td><td>0</td><td>0.0</td><td>Spatially varying (pysm d1s1)</td></tr><tr><td>3</td><td>Random for each band</td><td>0.0</td><td>Spatially varying (pysm d1s1)</td></tr><tr><td>4</td><td>Random for each band</td><td>0.3</td><td>Spatially varying (pysm d1s1)</td></tr></tbody></table>

by sets of non-linear parameters $\beta_d$, $T_d$ (parameterizing a modified blackbody) and $\beta_s$ (parameterizing the synchrotron's power law). The analytical parameterization of the SEDs in CMB units are given by

$$ 
\begin{align} 
\mathbf{A}_{\text{CMB}}^{X}(\nu, \nu_0) &= 1, \nonumber \\ 
\mathbf{A}_{\text{dust}}^{X}(\nu, \nu_0, \beta_d, T_d) &= \frac{1}{u(\nu)} \left(\frac{\nu}{\nu_0}\right)^{\beta_d - 2} \frac{B_{\nu}(T_d)}{B_{\nu_0}(T_d)}, \tag{2.3} \\ 
\mathbf{A}_{\text{sync}}^{X}(\nu, \nu_0, \beta_s) &= \frac{1}{u(\nu)} \left(\frac{\nu}{\nu_0}\right)^{\beta_s}, \nonumber 
\end{align} 
$$

where $X \in \{I, Q, U\}$, $B_\nu$ is the Planck blackbody function, and $u(\nu) = x^2 e^x / (e^x - 1)$, with $x = h\nu/(k_B T_{\text{CMB}})$, is the unit conversion factor transforming from antenna temperature to thermodynamic units. For a given frequency band $f$ and a given Stokes parameter, $\mathbf{A}_{\text{CMB}}$, $\mathbf{A}_{\text{dust}}$, and $\mathbf{A}_{\text{sync}}$ are simply three scalars. The quantity $s_x$ from eq. (2.1) are the component maps for the CMB, dust, and synchrotron, defined at the reference frequency $\nu_0 = 150$ GHz, while $n$ is the additional Gaussian white noise.

In the following, we describe these simulations in detail, with properties summarized in table 1.

* **Phase 1:** these are made of CMB maps generated from the Planck2018 best fit [61], combined with white noise based on the LiteBIRD latest instrument model [14], as well as foreground emission based on the pysm<fnref n="5" /> d0s0 model [62]. In this case, dust and synchrotron are considered “simple” because their SEDs are spatially constant.

* **Phase 2:** these are similar to Phase 1 but this time the foreground emission is based on the d1s1 model from pysm: the dust and synchrotron SEDs are spatially varying, as estimated from the analysis of Planck's total intensity data with Commander [63].

* **Phase 3:** these are built on top of Phase 2, with the inclusion of randomly distributed instrumental miscalibrations following ref. [60]. Each of the 100 simulations corresponds to independent realisations of the $\alpha_f$ parameter that creates the rotation matrix $R^{\text{inst}}$ in eq. (2.1). Similarly to ref. [64], the $\alpha_f$ values are drawn from independent uniform distributions centered at zero and with a width $\alpha_f \in [-\sigma_{\alpha_f}, \sigma_{\alpha_f}]$, where the $\sigma_{\alpha_f}$ for each frequency band is taken from “Case 2.3” summarized in table 3 of [60]

<footnote n="5">⁵https://github.com/galsci/pysm</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Average offset across frequency bands computed from the 100 simulations, before ((1/22) × Σ<sub>f</sub> α<sub>f</sub>) and after (Σ<sub>f</sub> W<sub>f</sub>α<sub>f</sub>, with the W<sub>f</sub> weights defined in section 3.1) component separation.</figcaption></figure>

<table><caption>**Table 2.** α<sub>f</sub> requirements for all frequency bands in the Low-, Mid-, and High-Frequency Telescopes (LFT, MFT, HFT) of *LiteBIRD*. These requirements on *LiteBIRD*'s absolute angle calibration per band ensure an unbiased measurement of the tensor-to-scalar ratio when grouping detectors per frequency and assuming that all frequencies are uncorrelated (Case 2.3 from table 3 of ref. [60]).</caption><thead><tr><th>Instrument</th><th>LFT</th><th>LFT</th><th>LFT</th><th>LFT</th><th>LFT</th><th>LFT</th><th>LFT</th><th>LFT</th><th>LFT</th><th>LFT</th></tr></thead><tbody><tr><td>Frequency [GHz]</td><td>40</td><td>50</td><td>60</td><td>68a</td><td>68b</td><td>78a</td><td>78b</td><td>89a</td><td>89b</td><td>100</td></tr><tr><td>σ<sub>α<sub>f</sub></sub> [arcmin]</td><td>147.8</td><td>118.2</td><td>47.9</td><td>41.8</td><td>106.5</td><td>25.5</td><td>38.6</td><td>15.9</td><td>87.4</td><td>11.3</td></tr></tbody></table><table><thead><tr><th>Instrument</th><th>LFT</th><th>MFT</th><th>MFT</th><th>MFT</th><th>MFT</th><th>MFT</th><th>HFT</th><th>HFT</th><th>HFT</th><th>HFT</th><th>HFT</th></tr></thead><tbody><tr><td>Frequency [GHz]</td><td>140</td><td>100</td><td>119</td><td>140</td><td>166</td><td>195</td><td>195</td><td>235</td><td>280</td><td>337</td><td>402</td></tr><tr><td>σ<sub>α<sub>f</sub></sub> [arcmin]</td><td>5.6</td><td>7.6</td><td>3.4</td><td>4.3</td><td>3.3</td><td>5.3</td><td>11.6</td><td>12.3</td><td>20.1</td><td>50.9</td><td>237.6</td></tr></tbody></table>

that we reproduce in table 2. For each simulation, the average miscalibration across frequencies, ⟨α<sub>f</sub>⟩<sub>f</sub>, gives a non-zero global offset to the whole instrument by a statistical fluctuation roughly within the limits considered in table 2 of ref. [60]. Figure 1 shows the distribution of average offsets before ⟨α<sub>f</sub>⟩<sub>f</sub> ≡ (1/22) × Σ<sub>f</sub> α<sub>f</sub>) and after component separation (⟨α<sub>f</sub>⟩<sub>f</sub> ≡ Σ<sub>f</sub> W<sub>f</sub>α<sub>f</sub>, where the weights W<sub>f</sub> will be defined in section 3.1) across the 100 simulations. The smaller spread of the weighted ⟨α<sub>f</sub>⟩<sub>f</sub> shows how the combination of frequency bands during component separation reduces the impact of relative miscalibrations.

* **Phase 4:** these are built on top of Phase 3 with the inclusion of an isotropic β = 0.3°, an amplitude that is representative of what was found in refs. [49–52].

We verified that, at the angular scales relevant for our analysis, *pysm's* d0 and d1 models contain a small non-null *EB* signal that is roughly compatible with *Planck's* dust
measurements within the statistical uncertainty of *Planck* polarisation data. Although enough
to introduce a noticeable bias in our analysis, such *EB* correlation is not purposely built into
d0 and d1 to reproduce the physical properties of dust emission. Hence, our simulations likely
underestimate the true foreground complexity that can arise from the coupling of magnetic
misalignments [53, 54, 65] and the spatial variation of the emission properties [66–69]. We
leave the study of more complex dust models (e.g., [70–73]) for future work.

## 3 Techniques employed

We employ five different pipelines to estimate $\beta$ in the presence of foreground contamination and systematic effects, as described in section 2. These pipelines are categorized into two groups. The first category, detailed in section 3.1, includes two pipelines that neither characterise nor control the effect of the instrumental $\alpha_f$. The second category, discussed in section 3.2, comprises three pipelines: one that implements the MK technique and two that fit $\alpha_f$ while performing component separation.

In order to improve readability, before delving into the details of each technique, we provide a concise description of the pipelines used to deliver the forecasts on CB. This is done through figure 2, which highlights the main features of each pipeline. In the left part of this figure, the inputs considered (shown in light grey) are divided into four different sets of simulations (i.e., Phases). Each of these Phases has been analysed by five pipelines, whose names are displayed in white boxes in figure 2. The pipelines in the blue area do not account for $\alpha_f$, while those in the orange area do. The former are applied to the CMB solutions provided by FGBuster, whereas the latter either extend their component separation to handle $\alpha_f$ or have no component separation at all since the method is applied directly to the frequency maps. In figure 2, below each pipeline name, the section where the pipeline is described in detail is indicated, while a brief description of the pipeline is provided next to it.

### 3.1 Pipelines without control of instrumental polarisation angles

These pipelines are both based on the CMB solution provided by the considered component-separation layer, i.e. FGBuster,<fnref n="6" /> a maximum-likelihood parametric map-based tool [74]. In order to increase the signal-to-noise ratio of the CB effect, we need to access high CMB multipoles. This is possible by neglecting the frequency channels that do not possess sufficiently high angular resolution. However, this choice has to be balanced against our ability to remove the foreground emission. The low-frequency channels have the lowest angular resolution, but trace the synchrotron emission the best. Hence, as a trade-off, we consider only medium and high-frequency channels with a full width at half maximum (FWHM) smaller than 30 arcmin. In figure 3, we display the signal-to-noise ratio we expect from this choice. In figure 4, we show that the bias induced by the synchrotron emission is subdominant in our results.

We proceed in two steps.

1. As in ref. [14], we degrade the resolution of the input, frequency-dependent maps, eq. (2.1), to a common 80 arcmin resolution. This procedure allows an unbiased

<footnote n="6">⁶https://github.com/fgbuster/fgbuster</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Overview of the pipelines employed to estimate CB. On the left, “wn” stands for white noise, d0s0 and d1s1 represent the foreground models, and α indicates the instrumental polarisation angles (see also section 2). On the right, the pipelines’ names are displayed in white boxes accompanied by a brief description. We indicate below each pipeline the section where a full description can be found. Pipelines in the blue block do not account for αf, while those in the orange block do.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Forecasted cumulative signal-to-noise on β from EB correlation (see eq. (1.6)) taking into account cosmic variance, noise, and foreground residuals from Phase 1 as estimated from FGBuster. The sky fraction is assumed to be f_sky = 60%, and l_min = 50.</figcaption></figure>

determination of the spectral indices, which would otherwise suffer from frequency-
dependent resolution. We then optimise the so-called spectral likelihood to estimate
the non-linear parameters of the synchrotron and dust SEDs, βs, βd and Td in eq. (2.1).

2. We build a weighting operator $\mathbf{W} = (\mathbf{A}^T\mathbf{N}^{-1}\mathbf{A}) \mathbf{A}^T\mathbf{N}^{-1}$, estimated using the noise covariance and the SEDs estimated in step (1), which we use to estimate the foreground-
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Significance of the averaged bias on the reconstructed angular power spectra. The low significance of fluctuations confirms that the choices made (such as neglecting the synchrotron contribution) are valid.</figcaption></figure>

cleaned Q and U CMB maps $\tilde{s}_{CMB}$ as well as foreground-cleaned CMB spectra. Using the notation from eq. (2.1), the estimated sky signal can be written as

$$ \tilde{s}_{kc}^{\text{fgs-clean}} = \sum_{f} \mathbf{W}_{cf} d_{kf}, \qquad (3.1) $$

where $d_{kf}$ corresponds to the input, noisy, and foreground-contaminated frequency maps (not convolved to a common resolution), $k$ denotes the Stokes parameter and $c$ is the index of the sky components (CMB, dust, synchrotron). Similarly, the component-separated spectra can be estimated as

$$ \tilde{C}_{\ell}^{cc', \text{fgs-clean}} = \sum_{f, f'} \mathbf{W}_{cf}^{T} C_{\ell}^{ff'} \mathbf{W}_{c'f'}, \qquad (3.2) $$

where $C_{\ell}^{ff'}$ are the frequency cross-power spectra estimated from the input maps. Note that we apply the $f_{sky} = 60\%$ Galactic mask derived by the Planck collaboration<fnref n="7" /> to the foreground-clean maps, eq. (3.1), prior to the estimation of $C_{\ell}^{ff'}$. Finally, we only consider the CMB solutions to eqs. (3.1) and (3.2), i.e. $c = c' = \text{CMB}$.

The weights constructed in this approach are also applied to foreground-only, noise-only, and foregrounds+noise maps to estimate the noise and statistical-foreground-residuals bias [75] present in the reconstructed, foreground-cleaned CMB map. In the case of real data sets, we would rely on similar foreground + noise simulations or use data splits. As noted in ref. [76], statistical foreground-residuals do not depend on the actual amplitude of foregrounds, but rather on the number of degrees of freedom considered in the component separation, making such an approach robust in the case of complex foregrounds.

The CMB solution $\tilde{s}^{\text{fgs-clean}}$ coming out from FGBuster is considered as the input of two estimators, one is harmonic-based, using the so-called “D-estimator” defined on the observed polarised spectra [77, 78] (section 3.1.1), and the other is pixel-based and consists of the

<footnote n="7">⁷HFI_Mask_GalPlane-apo0_2048_R2.00.fits from https://pla.esac.esa.int/#maps.</footnote>
stacking of Q and U maps at the location of the T and E peaks [35, 79, 80] (section 3.1.2). Because of their different domains, in principle, these two independent estimators respond to residual systematic effects differently. Hence, having these two compatible estimates is a benefit from the point of view of the robustness of the analysis. However, in these pipelines, the estimate of β is degenerate with α$_{f}$ . This means that the total error budget has to be complemented with additional (and independent) information from the instrument calibration. Reversing the argument, forecasts based on these pipelines provide a reference to assess requirements on the knowledge of α$_{f}$ . This can be achieved by requiring that the uncertainty on α$_{f}$ be subdominant with respect to the statistical uncertainty of β .

3.1.1 Harmonic-based, D -estimator

To build the D -estimator, we need a model for the observed, foreground-cleaned CMB (see eq. (3.2)). From eqs. (1.2) and (1.3) and from eqs. (1.4) to (1.6) we can build harmonic-based estimators based on the observed T B and E B spectra, respectively. Since the former has a signal-to-noise ratio smaller than the latter [35, 78], throughout this work we employ only D EB,obs ℓ , which is defined as

$$
D_{\ell}^{EB,\text{obs}}(\hat{\beta}) \equiv C_{\ell}^{EB,\text{obs}} \cos(4\hat{\beta}) - \frac{1}{2}(C_{\ell}^{EE,\text{obs}} - C_{\ell}^{BB,\text{obs}}) \sin(4\hat{\beta}), \quad (3.3)
$$

where ˆ β is the estimate for the birefringence angle. The main property of D EB,obs ℓ ( ˆ β ) is that its expectation value goes to zero when ˆ β = β , i.e.,

$$
\langle D_{\ell}^{EB,\text{obs}}(\hat{\beta} = \beta) \rangle = 0, \quad (3.4)
$$

see e.g. ref. [78] for further details. Therefore, this suggests that it is enough to look for the zeros of D EB,obs ℓ ( ˆ β ) to find β [35, 77, 78]. This is typically done through standard χ $^{2}$-minimisation techniques.

In the following, we will consider the observed spectra to be the noise-debiased, estimated, foreground-cleaned CMB spectrum:

$$
C_{\ell}^{XY, \text{obs}} \equiv \tilde{C}_{\ell}^{XY, \text{fgs-clean}} - \tilde{N}_{\ell}^{XY}, \qquad (3.5)
$$

where ˜ N$_{ℓ}$ is the estimated noise bias combined with statistical foreground residuals, and ˜ C XY,fgs-clean ℓ the CMB × CMB component-separated spectra defined in eq. (3.2), with X, Y ∈ { E, B } . In addition, we perform a complementary analysis without any noise and foreground residual debiasing, i.e., assuming C XY,obs ℓ = ˜ C XY,fgs-clean ℓ .

3.1.2 Pixel-based estimator

Our second estimator, which also does not attempt to remove any miscalibration angle, is a pixel-based approach. This method consists of the stacking of Q and U maps at the peak location of the E -mode and temperature maps [22, 35, 79–81]. At each of the T and E extrema, we calculate the transformed Stokes parameters,

$$
\begin{align}
Q_r(\theta) &= -Q(\theta) \cos(2\phi) - U(\theta) \sin(2\phi), \nonumber \\
U_r(\theta) &= Q(\theta) \sin(2\phi) - U(\theta) \cos(2\phi), \tag{3.6}
\end{align}
$$
where $\phi$ is the angle from local east, with north pointing towards the Galactic north pole (anticlockwise), and $\theta$ is a distance radially outward from the center of the extrema, introduced in ref. [82]. These can be interpreted as a local E and B transformation, where the $Q_r(\theta)$ measures the tangential and radial components of the polarisation at the radial distances $\theta$ from the peaks, and $U_r(\theta)$ the polarisation at 45° to the tangential and radial directions. The expected profiles of the transformed Stokes parameters are given for T extrema by

$$ \begin{aligned} \langle Q_r^T \rangle (\theta) &= -\int \frac{\ell d\ell}{2\pi} W_{\ell}^T W_{\ell}^P (\bar{b}_{\nu} + \bar{b}_{\zeta} \ell^2) C_{\ell}^{TE,\text{obs}} J_2(\ell\theta), \\ \langle U_r^T \rangle (\theta) &= -\int \frac{\ell d\ell}{2\pi} W_{\ell}^T W_{\ell}^P (\bar{b}_{\nu} + \bar{b}_{\zeta} \ell^2) C_{\ell}^{TB,\text{obs}} J_2(\ell\theta), \end{aligned} \quad (3.7) $$

and for E extrema by

$$ \begin{aligned} \langle Q_r^E \rangle (\theta) &= -\int \frac{\ell d\ell}{2\pi} W_{\ell}^E W_{\ell}^P (\bar{b}_{\nu} + \bar{b}_{\zeta} \ell^2) C_{\ell}^{EE,\text{obs}} J_2(\ell\theta), \\ \langle U_r^E \rangle (\theta) &= -\int \frac{\ell d\ell}{2\pi} W_{\ell}^E W_{\ell}^P (\bar{b}_{\nu} + \bar{b}_{\zeta} \ell^2) C_{\ell}^{EB,\text{obs}} J_2(\ell\theta), \end{aligned} \quad (3.8) $$

where $W_{\ell}^X$ with $X = P, E$ and $T$ are the beam and window pixel functions for the P polarisation (observed Q and U maps), E, and T maps respectively, $J_2(\ell\theta)$ are the second-order Bessel functions of the first kind, and $\bar{b}_{\nu}$ and $\bar{b}_{\zeta}$ are the scale-dependent bias parameters as described in ref. [79]. The E-mode map is produced from the observed Q and U maps at the same resolution as the T, Q, and U maps. While this does not have to be the case, we found no improvements with additional smoothing of the E-mode map, so we maintain the input resolution. In the case of no parity-violating terms, we expect $\langle U_r^X \rangle(\theta)$ to be zero (as there is no $C_{\ell}^{EB}$ or $C_{\ell}^{TB}$ predicted by $\Lambda$CDM). However, with CB modifying the power spectra as given in eqs. (1.3) and (1.6), we instead have,

$$ \begin{aligned} \langle U_r^T \rangle (\theta) &= -\sin(2\beta) \int \frac{\ell d\ell}{2\pi} W_{\ell}^T W_{\ell}^P (\bar{b}_{\nu} + \bar{b}_{\zeta} \ell^2) C_{\ell}^{TE} J_2(\ell\theta), \\ \langle U_r^E \rangle (\theta) &= -\frac{1}{2} \sin(4\beta) \int \frac{\ell d\ell}{2\pi} W_{\ell}^E W_{\ell}^P (\bar{b}_{\nu} + \bar{b}_{\zeta} \ell^2) (C_{\ell}^{EE} - C_{\ell}^{BB}) J_2(\ell\theta), \end{aligned} \quad (3.9) $$

which is different from zero for non-null $\beta$. It is then simple to perform a least-squares analysis on the resulting stacked profiles compared with the theory to obtain the final $\beta$. See refs. [22, 35, 81] for similar analyses and details.

## 3.2 Pipelines with control of instrumental polarisation angles

### 3.2.1 Template-based extension of the Minami-Komatsu estimator

This pipeline is based on the MK estimator introduced in refs. [46–48] and applied to data in refs. [49–52]. As we mentioned in the introduction, thanks to the use of the information contained in the foreground emission it is possible to estimate at the same time $\beta$ and $\alpha_f$ corresponding to each frequency channel $\nu_f$ (with $f = 1, \dots, 22$ in the case of LiteBIRD). This technique is robust against several instrumental systematics [55], including the miscalibration of $\alpha_f$ and other systematic effects that produce spurious EB correlations (e.g., intensity-to-polarisation leakage, beam leakage, or cross-polarisation effects [83–85]). However, it is
sensitive to the EB correlation of Galactic foreground emission and its modelling [50, 53–55, 70]. Hence, compared to the pipelines from section 3.1, this approach has the advantage of being independent of additional knowledge of the instrument calibration, at the price of needing a precise model for the EB correlation of the polarised foreground emission.

The MK pipeline is applied directly on frequency maps rather than the CMB solution provided by the component-separation layer to exploit the information contained in foreground emission. Taking foregrounds into account, the EB spectrum obtained from the cross-correlation of the *i* and *j* frequency channels is, within the small-angle approximation ($|α_f| \lesssim 5°$ and $|β| \lesssim 5°$),

$$C_{\ell}^{E_i B_j, \text{obs}} \approx 2\alpha_j C_{\ell}^{E_i E_j, \text{obs}} - 2\alpha_i C_{\ell}^{B_i B_j, \text{obs}} + 2\beta (C_{\ell}^{E_i E_j, \text{CMB}} - C_{\ell}^{B_i B_j, \text{CMB}}) \\ + \mathcal{A}^{\text{sync}} C_{\ell}^{E_i B_j, \text{sync}} + \mathcal{A}^{\text{dust}} C_{\ell}^{E_i B_j, \text{dust}} \\ + \mathcal{A}^{\text{sync}\times\text{dust}} (C_{\ell}^{E_i^{\text{sync}} B_j^{\text{dust}}} + C_{\ell}^{E_i^{\text{dust}} B_j^{\text{sync}}}), \quad (3.10)$$

where the superscripts indicate whether the spectra come from the correlation of the *E* and *B* modes of the observed frequency maps, the synchrotron and dust models, or the CMB spectra in the case of null CB, multiplied by the beam and pixel window function of each channel. The $\mathcal{A}^x$ are multiplicative factors added *ad hoc* to control the amplitude of the synchrotron ($\mathcal{A}^{\text{sync}}$) and dust ($\mathcal{A}^{\text{dust}}$) contributions as well as their cross-correlation ($\mathcal{A}^{\text{sync}\times\text{dust}}$). This constitutes an extension of the methodology compared to previous works [50, 55] where neither the synchrotron emission nor its cross-correlation with dust were considered. Although eq. (3.10) is a general expression considering a potential EB contribution from both synchrotron and dust, no significant synchrotron EB correlation has been found yet [86, 87] and we only have moderate evidence of non-null dust EB correlation induced by the misalignment between dust filaments and the Galactic magnetic field [53, 54, 65].

Following refs. [50, 55], we leave $\mathcal{A}^x$ as a free amplitude parameter and obtain the synchrotron and dust EB spectrum from the templates produced in component-separation analyses that fit parametric models to CMB data (e.g., Commander [88], FGBuster [74], or B-SeCRET [89]). In particular, we use the d0s0 and d1s1 foreground models provided by pysm [62, 90]. We take $\mathcal{A}^x$ as a single overall amplitude and use the SEDs provided by the parametric models to scale the foreground template to the target frequencies. See refs. [46, 47, 50, 51, 54, 70] for alternative ways to model the foreground EB correlation.

From eq. (3.10), we build a Gaussian likelihood to simultaneously fit ($\mathcal{A}^x, \beta, \alpha_f$) and find the maximum-likelihood solution for all parameters following the semi-analytical iterative algorithm presented in refs. [55, 91]. To emulate real analyses [49–52], we consider a minimal mask covering extragalactic point sources and the regions of brightest CO emission that leaves a sky fraction of $f_{\text{sky}} = 93.2\%$ free for the analysis after applying a 0.5° smoothing apodisation. For brevity, we limit our study to only one mask, since a good correction of the dust-EB bias is possible independently of the choice of mask when the foreground model that generated the simulations is also used as the dust template in the MK estimator [55]. In section 4.2.1, we also explore how our imperfect knowledge of foreground emission impacts the efficiency of the MK estimator by using templates of a simpler sky model to describe
a more complex sky. We use NaMaster<fnref n="8" /> [92] to calculate the pseudo-$C_ℓ$ of the masked sky and bin the power spectra into $N_b = 21$ uniform bins of $\Delta ℓ = 20$ in the range $ℓ \in [21, 440]$. In the presence of an unexpectedly high $1/f$ correlated noise contribution or an unreliable dust EB model, the minimum multipole can be raised to $ℓ_{min} \approx 100$ values with minimal impact on the results presented in section 4.2.1. Extending the analysis to higher multipoles does not report an advantage as the MK-estimator signal-to-noise on $\beta$ already saturates around $ℓ_{max} \approx 400$ for LiteBIRD's instrumental configuration. See appendix A for more details on the definition of the estimator.

### 3.2.2 CAB-SECRET and $\beta$CAB-SECRET

CAB-SECRET is shorthand for Calibration Angles B-SECRET, representing an extension of the standard B-SECRET method by incorporating $\alpha_f$ into the component-separation process. The component-separation analysis is performed using B-SECRET (Bayesian-Separation of Components and Residuals Estimate Tool) [89], an alternative parametric component-separation approach to FGBuster. B-SECRET utilizes an affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) to sample the posterior distribution, from which we can recover the marginalized parameters' posterior distribution. By integrating $\alpha_f$ into the component-separation process, they are marginalized over, and their uncertainties are propagated into the CMB maps, and hence into the estimation of $\beta$. Instead of sampling the entire posterior distribution simultaneously, we partition the parameter space into three subspaces: amplitudes ($A$); spectral parameters ($B$); and rotation angles ($C$). We then sample their conditional distributions. This iterative process continues until convergence is reached.

B-SECRET is a Bayesian method hence prior information is applied. Specifically, multivariate Gaussian priors are employed for $\alpha_f$ ($N(\alpha_f, C_{\alpha_f})$). These priors are derived from multi-frequency data by nulling the EB power spectra, as described in refs. [64, 91]. However, this method assumes $\beta = 0$. Consequently, if $\beta$ is non-zero, the recovered $\alpha_f$ become $\alpha_f + \beta$. In such cases, $\alpha_f$ estimates from the template-based MK estimator, as detailed in section 3.2.1, are utilized as prior information. The former set of priors is denoted as $\Pi(\alpha_f, \beta = 0)$, while the latter is denoted as $\Pi(\alpha_f, \beta)$. Alternatively, other prior information, such as instrument calibration before launch, external calibration with Tau A, or an artificial calibration source, could potentially be utilized.

Prior to conducting the component-separation analysis, the computational load for fitting the maps at $N_{side} = 512$ is reduced by downgrading them to $N_{side} = 64$. Additionally, all maps are subjected to convolution with the same beam as follows.

(i) The original maps at $N_{side} = 512$ are transformed into spherical harmonics coefficients.

(ii) Pixel window function and beam deconvolution is applied to each frequency channel based on the beam's FWHM, as specified in table 3 of ref. [14].

(iii) The spherical harmonic coefficients of each channel are convolved with a Gaussian beam with FWHM = 2.2° and the corresponding pixel window function, and then transformed into $N_{side} = 64$ maps. The FWHM corresponds to 2.4 times the effective pixel size at $N_{side} = 64$ to account for the pixel window function effects.

<footnote n="8">⁸https://github.com/LSSTDESC/NaMaster</footnote>
Additionally, we apply the same Planck $f_{\text{sky}} = 60\%$ Galactic-plane mask discussed in section 3.1 to mitigate the influence of large residuals from foregrounds on other parameters.

We consider two cases of this pipeline, **CAB-SECRET**, where $\alpha_f$ are included in the parametric model and $\beta$ is estimated from the cleaned CMB map, and **$\beta$ CAB-SECRET**, where $\beta$ is included as another parameter in $C$.

**CAB-SECRET.** CAB-SECRET includes only the $\alpha_f$ as model parameters ($C = \{\alpha_f, \forall f \in N_\nu\}$) as in ref. [91]. Hence, the signal is modelled as follows for a given pixel $p$:

$$ \begin{pmatrix} Q^{\text{obs}}(f, \alpha_f, \gamma) \\ U^{\text{obs}}(f, \alpha_f, \gamma) \end{pmatrix}_p = \begin{pmatrix} \cos(2\alpha_f) & -\sin(2\alpha_f) \\ \sin(2\alpha_f) & \cos(2\alpha_f) \end{pmatrix} \begin{pmatrix} Q(f, \gamma) \\ U(f, \gamma) \end{pmatrix}_p = \mathbf{R}(\alpha_f) \begin{pmatrix} Q(f, \gamma) \\ U(f, \gamma) \end{pmatrix}_p, \quad (3.11) $$

where $\mathbf{R}(\alpha_f)$ is the rotation matrix that accounts for the $Q$ and $U$ mixing due to non-perfect calibration, $\gamma = A \cup B$, and

$$ \begin{pmatrix} Q(f, \gamma) \\ U(f, \gamma) \end{pmatrix}_p = \begin{pmatrix} a_{\text{CMB}}^Q \\ a_{\text{CMB}}^U \end{pmatrix}_p + \begin{pmatrix} a_s^Q \\ a_s^U \end{pmatrix}_p \frac{1}{u(f)} \left(\frac{f}{f_s}\right)^{\beta_s} + \begin{pmatrix} a_d^Q \\ a_d^U \end{pmatrix}_p \frac{1}{u(f)} \left(\frac{f}{f_d}\right)^{\beta_d - 2} \frac{B_f(T_d)}{B_{f_d}(T_d)}, \quad (3.12) $$

is the parametric model used to describe the astrophysical components considered in the simulations, taking $f_s = 30$ GHz and $f_d = 353$ GHz as pivot frequencies for the synchrotron and dust emission.

The amplitude parameters

$$ A = \bigcup_{p \in N_{\text{pix}}} \{a_{\text{CMB}}^Q, a_{\text{CMB}}^U, a_s^Q, a_s^U, a_d^Q, a_d^U\}_p \quad (3.13) $$

vary from pixel to pixel, while the spectral parameters

$$ B = \bigcup_{c \in N_c} \{\beta_s, \beta_d, T_d\}_c \quad (3.14) $$

are considered uniform within specific sky regions (i.e., clusters of pixels). Additionally, individual Gaussian priors are applied on $B$: $\beta_s \sim N(-3.1, 0.3)$, $\beta_d \sim N(1.56, 0.1)$, and $T_d \sim N(21, 3)$ [89, 91]. Reducing the number of parameters by clustering the $B$ yields significantly smaller statistical uncertainties for all parameters, including CB, although it introduces biases in $\beta$ when the assumption of constant $B$ is invalid. Different clustering strategies are investigated in section 4.2.2 using Phase 2 simulations where spatial variations of the spectral parameters are present.

From the cleaned CMB maps, we determine CB as the residual rotation angle in the CMB, following the approach outlined in ref. [64]. This involves finding the $\beta$ value that minimises the harmonic estimator derived from eq. (3.3). Initially, we calculate the angular power spectrum and their covariance matrix using NaMaster. The power spectra are then binned uniformly with a bin width of $\Delta\ell = 10$ within the range $\ell \in [12, 191]$, employing the same weighting scheme as ref. [64]. We choose to use all the available bins at $N_{\text{side}} = 64$, except the first bin due to the high uncertainties recovered with a pseudo-$C_\ell$ method. Additionally,
we apply apodisation to the mask using a C2 kernel with an apodisation scale of 22°, as suggested in ref. [93].

The CB angle is obtained by minimising the following log-Gaussian distribution:

$$-2 \log \mathcal{L} = D^{\text{EB,obs}}(\beta)\mathbf{M}^{-1}D^{\text{EB,obs}}(\beta) + \log |\mathbf{M}| + N_b \log (2\pi), \quad (3.15)$$

where $D^{\text{EB,obs}}(\beta)$ is the vector of harmonic-based estimators in eq. (3.3) evaluated for each multipole bin $b$, $N_b$ is the number of multipole bins, $|...|$ is the determinant operator, and

$$\mathbf{M} = \cos^2(4\beta)\mathbf{\Sigma}^{\text{EB,EB}} + \cos(4\beta)\sin(4\beta)(\mathbf{\Sigma}^{\text{EB,BB}} - \mathbf{\Sigma}^{\text{EB,EE}}) \\ + \frac{1}{4}\sin^2(4\beta)(\mathbf{\Sigma}^{\text{EE,EE}} + \mathbf{\Sigma}^{\text{BB,BB}} - 2\mathbf{\Sigma}^{\text{EE,BB}}), \quad (3.16)$$

where $\mathbf{\Sigma}$ represents the covariance matrix of the cross-power spectra, i.e.,

$$\Sigma_{bb'}^{XY,X'Y'} \equiv \langle \Delta C_b^{XY} \Delta C_{b'}^{X'Y'} \rangle. \quad (3.17)$$

We employ emcee<fnref n="9" /> [94] to sample the likelihood in eq. (3.15). The median and standard deviation of the resulting chains provide estimates for $\beta$ and $\sigma_\beta$, respectively.

**βCAB-SECRET.** In this pipeline, CB is incorporated into the model as a parameter, i.e., $C = \{\beta\} \cup \{\alpha_f, \forall f \in N_\nu\}$. Thus, eqs. (3.11) and (3.12) lead to

$$\begin{pmatrix} Q^{\text{obs}}(f, \beta, \alpha_f, \gamma) \\ U^{\text{obs}}(f, \beta, \alpha_f, \gamma) \end{pmatrix}_p = \mathbf{R}(\beta + \alpha_f) \begin{pmatrix} a_{\text{CMB}}^Q \\ a_{\text{CMB}}^U \end{pmatrix}_p + \mathbf{R}(\alpha_f) \left[ \begin{pmatrix} Q(f, \gamma) \\ U(f, \gamma) \end{pmatrix}_p - \begin{pmatrix} a_{\text{CMB}}^Q \\ a_{\text{CMB}}^U \end{pmatrix}_p \right]. \quad (3.18)$$

In this scenario, the cleaned CMB map undergoes no additional rotation, since any rotation has been accounted for and removed during the component separation process. Consequently, $\beta$ is directly obtained from the component-separation step.

For βCAB-SECRET, we utilise the results from the template-based MK estimator as prior information, denoted as $\Pi(\alpha_f, \beta)$. It is important to note that this approach differs from the template-based MK pipeline in that we simultaneously allow for variations in the foreground model and CB. This simultaneous variation helps mitigate biases arising from potentially incorrect modelling of foreground templates, as demonstrated in section 4.2.2.

### 3.2.3 Generalised pixel-based parametric component separation, J23

Similarly to the previous section, the pipeline presented here includes the effect of $\alpha_f$ in its data model. This allows us to fit them in the component-separation step, and therefore correct their effects on the CMB map. This pipeline is a generalisation of FGBuster [74], which is the method used to obtain the cleaned CMB maps employed in sections 3.1.1 and 3.1.2. The method we use here is therefore pixel-based and parametric. The details of its generalisation for the inclusion of $\alpha_f$ and CB are described in ref. [95]. For conciseness, we will call this pipeline J23. The data model used to describe the frequency maps $\mathbf{d}_{hp}$ is as follows:

$$\mathbf{d}_{hp} = \mathbf{X}_{hh'}(\{\alpha_f\})\mathbf{A}_{h'o}(\{\beta_d, T_d, \beta_s\})\mathbf{s}_{op} + \mathbf{n}_{hp}. \quad (3.19)$$

Here $p$ describes pixels. Index $h$ spans both $\{Q, U\}$ Stokes parameters and frequency bands (the number of $h$ indices is then $n_h = n_{\text{Stokes}} n_{\text{frequency}}$). Similarly, $o$ covers $\{Q, U\}$ and sky components ($n_o = n_{\text{Stokes}} n_{\text{component}}$). $\mathbf{s}_{op}$ is the vector storing the $\{Q, U\}$ maps for

<footnote n="9">⁹https://github.com/dfm/emcee</footnote>
each considered sky component: {CMB, dust, synchrotron}. $\mathbf{X}_{hh'}$ is an $n_{\text{Stokes}} n_{\text{frequency}} \times n_{\text{Stokes}} n_{\text{frequency}}$ block diagonal matrix that encodes the effect of $\alpha_f$, each block being a rotation matrix $\mathbf{R}^{\text{inst}}(\alpha_f)$ acting on the $\{Q, U\}$ maps at frequency $f$. $\mathbf{A}_{ho}$ is the mixing matrix, encoding the frequency dependence of the different components, which is an $n_{\text{Stokes}} n_{\text{frequency}} \times n_{\text{Stokes}} n_{\text{component}}$ matrix. We assume: (i) the frequency scaling to be the same between the $Q$ and $U$ parameters of a given component; (ii) no intrinsic $Q-U$ mixing, making every other element of a column equal to 0; (iii) that the CMB scaling is known and that we are working in CMB units, thus the non-zero elements of the first two columns of $\mathbf{A}$ are all equal to 1; and (iv) that the dust frequency scaling is a standard modified blackbody law parameterised by the dust spectral index $\beta_d$ and the dust temperature $T_d$, and that the synchrotron is a power law with index $\beta_s$ (see eq. (2.3) for a description of those laws). Finally, $\mathbf{n}_{hp}$ is the white Gaussian noise present in the frequency maps, same as in eq. (2.1). For conciseness, we define $\mathbf{\Lambda}_{ho}(\beta_d, T_d, \beta_s, \{\alpha_f\}) \equiv \mathbf{X}_{hh'}(\{\alpha_f\})\mathbf{\Lambda}_{h'o}(\beta_d, T_d, \beta_s)$.

Once we have the data model we can proceed to the analysis. The pipeline goes through 5+1 steps that are detailed below.

**Step 0: pre-processing.** Contrary to eq. (2.1), the data model in eq. (3.19) does not explicitly include the effect of the beam. To deal with this, the frequency maps are first deconvolved according to their respective beam (table 3 of ref. [14]) and convolved to the same common Gaussian beam with FWHM of 80 arcmin used in ref. [14], as well as in the standard FGBuster approach in section 3.1. This is done so that we do not suffer from frequency dependence of the resolution. Furthermore, convolving to a smaller FWHM would boost the noise of channels with the poorest resolution. Extending the data model to include the beams at different frequencies is explored in ref. [96], and is a possible extension of this method, but beyond the scope of the work presented here. The maps are also downgraded to $N_{\text{side}} = 64$, deconvolved by the pixel window function corresponding to the old $N_{\text{side}} = 512$, and convolved with the new one. The data vector after pre-processing is denoted $\bar{\mathbf{d}}$. This operation will notably affect the white noise, which is then denoted $\bar{\mathbf{n}}$. We apply Planck's Galactic masks with either $f_{\text{sky}} = 60\%$ or $f_{\text{sky}} = 40\%$ depending on the case, as will be described in section 4.2.3.

**Step 1: component separation.** Unlike B-SeCRET, FGBuster and its generalisation use a ridge-likelihood, called spectral likelihood, $\mathcal{L}^{\text{spec}}$, where the only free parameters are $\{\beta_d, T_d, \beta_s, \{\alpha_f\}\}$ encoding $\mathbf{\Lambda}$. Here, we assume those parameters to be constant across the sky.

While it was already shown in ref. [95] that the J23 pipeline correctly retrieves the relative angles between frequencies, the degeneracy between the global $\alpha$ offset and $\beta$ remains. To lift it, additional information must be added to the likelihood. Here, it takes the form of Gaussian priors on the $\{\alpha_f\}$ parameters added to the likelihood. The generalized spectral likelihood that we use is then

$$ -2 \log(\mathcal{L}^{\text{spec}}(\{\beta_{\text{fg}}, \alpha_f\})) = + \text{constant} - \sum_{\text{p}} \text{tr} \left( \overline{\mathbf{N}}^{-1} \mathbf{\Lambda} (\mathbf{\Lambda}^{\top} \overline{\mathbf{N}}^{-1} \mathbf{\Lambda})^{-1} \mathbf{\Lambda}^{\top} \overline{\mathbf{N}}^{-1} \bar{\mathbf{d}}_{\text{p}} \bar{\mathbf{d}}_{\text{p}}^{\top} \right) + \sum_{\text{f}} \frac{(\alpha_{\text{f}} - \tilde{\alpha}_{\text{f}})^2}{\sigma_{\alpha_{\text{f}}}^2}, \quad (3.20) $$
with constant being the constant ensuring the positivity of the right-hand side. We have omitted the o and h indices for conciseness, and f spans the frequency bands. Although the pre-processing correlates the initial white noise, here we still assume a white and scale-invariant noise in our analysis. $\overline{N}$ is therefore a pixel-independent and diagonal $n_{\text{Stokes}} n_{\text{frequency}} \times n_{\text{Stokes}} n_{\text{frequency}}$ matrix. The diagonal elements of $\overline{N}$ are computed from $\langle \overline{n} \overline{n}^\top \rangle$, where we average over noise simulations and pixels. This approximate modelling of the noise properties does not seem to introduce a significant bias in our results, which, as will be discussed in section 4.2.3, are dominated by the approximations in foreground modelling.

The second term in eq. (3.20) corresponds to the Gaussian priors. We assume the priors on each $\alpha_f$ to be independent. The prior on the angle $\alpha_f$ is centred on $\tilde{\alpha}_f$ and has a variance of $\sigma_{\alpha_f}^2$. Our approach here is to assume priors coming from an independent ground or in-flight calibration. As such, we assume these priors match the current requirements of the LiteBIRD mission and take $\sigma_{\alpha_f}$ from table 2. The central values of the priors $\tilde{\alpha}_f$ are not chosen to be the true values of the input $\alpha_f$ but drawn from a normal distribution centred at zero with variance $\sigma_{\alpha_f}^2$.<fnref n="10" />

To estimate $\{\beta_d, T_d, \beta_s, \{\alpha_f\}\}$ we must minimise eq. (3.20). In practice, we initialise the minimisation by performing an MCMC sampling with emcee [94]. Once the MCMC converges, we use the mean of the chains to initialise the minimisation of eq. (3.20), which is performed using the L-BFGS-B method from scipy [97]. We denote $\hat{\beta}_d$, $\hat{T}_d$, $\hat{\beta}_s$, $\{\hat{\alpha}_f\}$ the best estimates of the generalized spectral likelihood that are obtained at the end of this process.

**Step 2: CMB map estimation.** With the best estimate of spectral parameters and $\alpha_f$ at hand, we can estimate the generalized mixing matrix $\hat{\Lambda}_{ho} \equiv \Lambda_{ho}(\hat{\beta}_d, \hat{T}_d, \hat{\beta}_s, \{\hat{\alpha}_f\})$. From there, we can build a weighting operator $\hat{W}_{oh}$ and apply it to the data similarly as in eq. (3.1):

$$ \hat{W}_{oh} = (\hat{\Lambda}_{oh'}^{\top} \overline{\mathbf{N}}_{h'h''}^{-1} \hat{\Lambda}_{h''o'}^{\top})^{-1} \hat{\Lambda}_{o'h'''}^{\top} \overline{\mathbf{N}}_{h'''h}^{-1}. \quad (3.21) $$

This operator can be used to retrieve all three components' maps, but we focus solely on its CMB part, $\hat{W}_{kh}^{\text{CMB}}$, where the k index represents the $\{Q, U\}$ Stokes parameters ($\hat{W}_{kh}^{\text{CMB}}$ shape is then $n_{\text{Stokes}} \times n_{\text{Stokes}} n_{\text{frequencies}}$). The $\{Q, U\}$ CMB maps after generalised component separation are then

$$ \hat{s}_{\text{CMB},kp}^{\text{fgs-clean}} = \hat{W}_{kh}^{\text{CMB}} \bar{d}_{hp}. \quad (3.22) $$

**Step 3: spectra estimation.** Once we have the cleaned CMB maps, we can estimate their angular power spectra using NaMaster. To do so, we use an apodised version of the mask used for component separation, based on a “smooth” apodization with a 5° size. The spectra are binned uniformly with $\Delta\ell = 2$ in the range $\ell \in [30, 125]$. In this step, the common beam of 80 arcmin is also corrected for. We denote the estimated CMB spectra $D_b$, with b representing the multipole at the centre of a given bin. $D_b$ is a $2 \times 2$ matrix with diagonal elements storing the EE and BB auto-spectra and the off-diagonal EB, assuming the EB cross-spectrum is equal to the BE one.

<footnote n="10">¹⁰In ref. [95], the calibration measurements were averaged over, leading to an additional factor of 2 in the variance and no bias in the prior.</footnote>
**Step 4: noise post-component separation.** To estimate the noise left in the angular power spectrum after all those steps, one needs to take into account the weighting operator, as well as the different beam convolution and deconvolution that have been applied to the white noise. After cancelling the effect of the common beam in the spectra estimation step, only the deconvolution is left. The noise power spectra in the CMB maps are then given by

$$C_b^{\text{noise}} = [\hat{\mathbf{W}}\mathbf{N}_b\hat{\mathbf{W}}^\text{T}]_{\text{CMB}}, \qquad (3.23)$$

where $\mathbf{N}_b$ is a $2n_{\text{frequency}} \times 2n_{\text{frequency}}$ matrix describing the noise power spectra for each frequency band. It is diagonal because we assume no $E-B$ correlation in the noise. Furthermore, we assume the noise has the same characteristics in $E$ and $B$. We can therefore describe the elements of $\mathbf{N}_b$ with respect to its bin and frequency $f$ as

$$\mathbf{N}_{b,f} = (w_f)^{-2} e^{b(b+1)\text{FWHM}_f^2/8 \log 2}. \qquad (3.24)$$

where $(w_f)^{-1}$ is the sensitivity at frequency $f$ in $\mu$K-rad and FWHM$_f$ the full width at half maximum of the beam of the same channel, in radians.

**Step 5: cosmological likelihood.** With the noise covariance post-component separation at hand, we can build a model covariance for the data:

$$\mathbf{C}_b^{\text{theory}}(\beta) = \mathbf{R}(\beta) \begin{pmatrix} C_b^{\text{EE, prim}} & 0 \\ 0 & C_b^{\text{BB, lens}} \end{pmatrix} \mathbf{R}^{-1}(\beta) + \mathbf{C}_b^{\text{noise}}. \qquad (3.25)$$

The first term describes the CMB power spectra. We used the Planck-2018 best-fit power spectra for the primordial CMB EE power spectrum, $C_b^{EE,prim}$, and for the lensed BB power spectrum, $C_b^{BB,lens}$. Contrary to ref. [95], the tensor-to-scalar ratio $r$ is not allowed to vary and is fixed to $r=0$. The BB power spectrum is hence limited to the lensing, and we fix $A_{\text{lens}} = 1$. We assume no intrinsic EB in the CMB, making the diagonal terms null. An EB signal can however appear via the effect of CB through the $\mathbf{R}(\beta)$ operators.

We can then compare our data $\mathbf{D}_b$ to this model covariance using the following Gaussian likelihood [98]:

$$-2 \log \mathcal{L}^{\text{cosmo}}(\beta) = \sum_b \Delta\ell f_{\text{sky}}(2b+1) \left[ \mathbf{C}_b^{\text{theory}}^{-1}(\beta) \mathbf{D}_b + \log(|\mathbf{C}_b^{\text{theory}}(\beta)|) \right]. \qquad (3.26)$$

We approximate the effective number of modes accessible after masking and binning with the use of the $\Delta\ell f_{\text{sky}}(2b+1)$ factor, therefore neglecting any mode coupling coming from the analysis of a partial sky. We can then proceed similarly as for the generalised spectral likelihood: first exploring the likelihood distribution with MCMC sampling, and then minimising eq. (3.26) starting from the mean of the converged MCMC chain. From there we obtain the final best estimate of $\beta$ for each simulation.
# 4 Analyses

## 4.1 Pipelines without control of instrumental polarisation angles

### 4.1.1 D-estimator analysis

As mentioned in section 3.1.1, CB is estimated looking for the $\hat{\beta}$ angle that minimises the following $\chi^2(\hat{\beta})$:

$$ \chi^2(\hat{\beta}) = \sum_{\ell, \ell'} D_{\ell}^{EB,obs}(\hat{\beta}) M_{\ell\ell'}^{-1} D_{\ell'}^{EB,obs}(\hat{\beta}), \quad (4.1) $$

where $M_{\ell\ell'} = \langle D_{\ell}^{EB,obs}(\beta) D_{\ell'}^{EB,obs}(\beta) \rangle$ is the covariance matrix of $D_{\ell}^{EB,obs}(\hat{\beta} = \beta)$, which is defined in eq. (3.3). It is possible to show that the matrix $M_{\ell\ell'}^{EB}$ reduces to the EB covariance matrix, i.e., without the CB effect included [78] when the noise of the E and B CMB fields are equivalent. Since this condition is always satisfied in all the cases we consider, we take $M_{\ell\ell'}$ to be the EB covariance. This covariance is built through simulations: we do that by performing an aggressive harmonic binning because of the low number of simulations available (details below). Moreover, for each of the simulations employed to sample the $\chi^2(\hat{\beta})$, the covariance $M_{\ell\ell'}$ is built from the remaining simulations (i.e., excluding the specific simulation considered in $D_{\ell}^{EB,obs}$). In the following, we discuss in detail the results obtained for each Phase.

* **Phase 1:** in figure 5, we show the spectra and the combination of the spectra needed to build $D_{\ell}^{EB,obs}$. Blue and orange curves are for the simulations in which CMB + noise + foregrounds and only noise + foregrounds are considered, respectively. Green is for the simulations in which the spectra have been debiased with noise + foregrounds. It is interesting to observe that the expected averages for these spectra are recovered in both the debiased and unbiased cases. Moreover, the foreground residuals show up at low multipoles, while the noise (coupled to the beam) presents itself at high multipoles. In between these regions, around $20 \le l \le 500$, the EB correlation is compatible with zero, which is expected in the absence of CB. For each simulation, we estimate $\beta$ and construct the histograms given in figure 6. These are obtained by building the covariance through simulations, setting a harmonic binning $\Delta l = 100$ in the range [100, 700], displayed in orange, or [50, 750], displayed in blue. In the same figure 6, we show histograms of $\beta$ obtained through spectra with and without noise + foregrounds debiasing. In particular, the latter are given with filled histograms, and the former are provided with unfilled histograms. In all the cases, the average of $\beta$ is found to be compatible with zero with an uncertainty of 0.011° at 68% C.L.. Moreover, approximating the covariance $M_{\ell\ell'}$ to be diagonal, we find the same statistical efficiency for all cases.

* **Phase 2:** considering the same setting adopted for Phase 1, we estimate $\beta$ for each simulation and construct the histograms given in figure 6. We find an uncertainty on $\beta$ of 0.011° at 68% C.L. for all the cases considered. Again similarly to Phase 1, when we approximate the covariance $M_{\ell\ell'}$ to be diagonal, we find the same statistical efficiency for all the cases considered.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Upper panel: EB spectra. Lower panel: $(EE - BB)/2$ spectra after component separation using FGBuster. In all the panels, Phase 1 simulations of CMB + noise + foregrounds are displayed in blue and simulations of only noise + foregrounds are shown in orange. The spectra debiased with foreground and noise are given in green. The shaded regions represent the standard deviation of the spectra across the 100 simulations.</figcaption></figure>

* **Phase 3:** following the same approach considered above, we construct the histograms given in figure 6. We find an uncertainty on $\beta$ of 0.046° at 68% C.L. for all the cases considered. Moreover, approximating the covariance $M_{\ell\ell'}$ to be diagonal, we find a general increase in the uncertainty on $\beta$ of a factor of 1.6 to 1.7. This seems to indicate that the inclusion of $\alpha_f$ different from zero (albeit compatible with requirements) provides an increase of the uncertainty of $\beta$ compared to Phases 1 and 2, and also non-negligible off-diagonal correlations, which, if not taken into account, make the uncertainty even larger.

* **Phase 4:** in figure 7, we show the spectra needed to build $D_{\ell}^{EB,obs}$. As in the previous cases, the averages of the recovered spectra are compatible with the fiducial spectra (which, in this case, are rotated with respect to the $\Lambda$CDM case by $\beta = 0.3^\circ$) for both the noise + foregrounds debiased and the non-debiased cases. Employing the same harmonic setting as above, we build the histograms shown in figure 6. When we consider the harmonic range $\ell \in [100, 700]$ we find $\beta = 0.296^\circ \pm 0.045^\circ$ and $\beta = 0.295^\circ \pm 0.045^\circ$ for the non-debiased and debiased cases, respectively. We hence retrieve the fiducial value of $\beta^{\text{fid}} = 0.3^\circ$. Similarly, for the harmonic range $\ell \in [50, 750]$, we obtain $\beta = 0.296^\circ \pm 0.047^\circ$ and $\beta = 0.297^\circ \pm 0.047^\circ$ for the non-debiased and noise + foregrounds debiased cases, respectively. Moreover, as for the Phase 3 case, we note that considering a diagonal matrix, we obtain an increase of the uncertainty on $\beta$ of a factor of 1.6 to 1.7.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Histograms of $\beta$ with and without debiasing, as obtained for all simulation phases. For each phase, we considered two harmonic ranges: [100, 700], shown in orange, and [50, 750], shown in blue. Shaded histograms represent the case without debiasing, while empty histograms indicate the case with debiasing.</figcaption></figure>

### 4.1.2 Pixel-based analysis

To determine $\beta$ using the pixel-based approach we begin by finding all the hot and cold extrema using the nearest-neighbour method on both *E* maps (produced from noisy *Q* and *U*), and the fiducial *T* maps used to generate the *LiteBIRD* simulations (clean with no noise). We chose a threshold of zero such that we only retain positive hot spots and negative cold spots. Including all peaks comes at the cost of more noise peaks adding no additional information to the analysis (increasing required computation time and slightly increasing the error bar), whereas including fewer peaks (with a more constraining threshold) comes at the risk of ignoring true peaks. Thresholds of zero have been used for previous analyses and we find that more constraining thresholds do not improve our results.

From eq. (3.9), we take $\theta$ out to 5°, where the profile converges towards zero. Higher angular sizes tested did not change the results. Peaks within 5° of the mask are discarded since they often lack enough pixels to both be identified as a peak and to correctly determine the value for $\beta$ (the sample of peaks near a masked region had poorer return of the fiducial $\beta$ and with much higher variance than other peak samples).
<figure><img src="image_10.png" /><figcaption>Figure 7. Same as figure 5 but for the Phase 4 simulations.</figcaption></figure>

<table><caption>**Table 3.** Stacking theory bias parameters for various data sets. Note that the bias parameters are very dependent on the map features (such as the beam, window function, and $N_{\text{side}}$), and the underlying power spectrum in the map (including the noise).</caption><thead><tr><th>Phase</th><th>N<sub>side</sub></th><th>Beam FWHM<br>[arcmin]</th><th>b<sub>ν</sub><sup>E</sup></th><th>b<sub>ζ</sub><sup>E</sup></th><th>b<sub>ν</sub><sup>T</sup></th><th>b<sub>ζ</sub><sup>T</sup></th></tr></thead><tbody><tr><td>1</td><td>64</td><td>80</td><td>0.226</td><td>2.85×10<sup>-4</sup></td><td>0.0022</td><td>1.47×10<sup>-6</sup></td></tr><tr><td>2</td><td>64</td><td>80</td><td>0.226</td><td>2.86×10<sup>-4</sup></td><td>0.0022</td><td>1.47×10<sup>-6</sup></td></tr><tr><td>3</td><td>128</td><td>40</td><td>0.157</td><td>2.67×10<sup>-5</sup></td><td>0.0018</td><td>4.06×10<sup>-7</sup></td></tr><tr><td>4</td><td>128</td><td>40</td><td>0.157</td><td>2.66×10<sup>-5</sup></td><td>0.0018</td><td>4.06×10<sup>-7</sup></td></tr></tbody></table>

For each peak, $m$, we obtain an estimate of the angle $\beta_m$, and a goodness of fit $\sigma^2_{\beta_m}$ using a standard least-squares analysis fit. First, we calculate the theoretical value $\tilde{U}_r^X(\theta_p)$, $X = \{T, E\}$, for all the pixels at $\theta < 5^\circ$ from the central peak, using eq. (3.9). The values for the bias parameters from the final maps used can be found in table 3, and the final theory profiles are shown in figure 8.

Next, we compute $\hat{U}_r(\theta_{\text{pix}})$ at the location of the $T$ or $E$ peak on the $Q$ and $U$ maps, using eq. (3.6). The peak values are calculated using

$$ \beta_m = \frac{1}{2} \frac{\sum_p \hat{U}_r(\theta_p) \tilde{U}_r^X(\theta_p)}{\sum_p \tilde{U}_r^X(\theta_p) \tilde{U}_r^X(\theta_p)}, \qquad (4.2) $$

$$ \sigma^2_{\beta_m} = \frac{1}{N-1} \frac{\sum_p (\hat{U}_r(\theta_p) - 2\beta_m \tilde{U}_r^X(\theta_p))^2}{\sum_p \tilde{U}_r^X(\theta_p) \tilde{U}_r^X(\theta_p)}, $$
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Theory profile for E and T peaks for different $N_{side}$, used in the pixel-based analysis. Phase 1 and 2 have $N_{side}$ = 64 and Phase 3 and 4 have $N_{side}$ = 128. Note that the theory profiles are shown here assuming $\beta = 0.3^{\circ}$, otherwise theory profiles for $\langle U_r^X \rangle(\theta)$ would be consistent with zero.</figcaption></figure>

where N is the total number of pixels in the sum, and $\sigma_{\beta_m}^2$ is the sum of the squares of the residuals and gives a goodness of fit to each peak calculation. Note that eq. (4.2) implicitly assumes that the noise covariance matrix around each peak is diagonal and independent of $\theta_m$. This approximation is enough to describe the statistics of the peaks’ surroundings in early simulation Phases but becomes insufficient for later Phases, and especially higher map resolutions, where component separation introduces more complicated residuals and pixel correlations. The construction of a more optimal estimator is left for future work.

We weight the $\beta_m$ using the residuals to obtain the final map result:

$$ \beta = \frac{\sum_m \beta_m / \sigma_{\beta_m}^2}{\sum_m 1 / \sigma_{\beta_m}^2}. \qquad (4.3) $$

These results are calculated on maps provided using the FGBuster method at $N_{side}$ = 64 with an FWHM beam of 80 arcminutes for all Phases, at $N_{side}$ = 128 with an FWHM beam of 40 arcminutes for Phase 3 and 4, and at $N_{side}$ = 512 with an FWHM beam of around 26 arcminutes (this is an approximation because the true beam for the highest resolution data is not a perfect Gaussian) also for Phases 3 and 4, as well as on the CAB-SECRET maps at $N_{side}$ = 64 with a beam of FWHM of 132 arcminutes for all Phases. The Galactic plane mask with $f_{sky}$ = 60% was applied to the data (as in pipeline CAB-SECRET and the harmonic approach).

Results for each of these maps are summarised in figure 9. The maps producing the best results came from the FGBuster method, for Phases 3 and 4 the $N_{side}$ = 128 results outperformed the $N_{side}$ = 512 results. Simple tests show that in maps with similar noise properties, increasing the resolution leads to better constraints on the results, however, in this case the $N_{side}$ = 512 maps have more complex noise properties. To combine the various frequency maps (with their distinct resolutions and beam functions) the method for producing the $N_{side}$ = 128 and the $N_{side}$ = 512 maps using FGBuster are very different.<fnref n="11" /> As a result, the

<footnote n="11">¹¹While the $N_{side}$ = 128 component maps were obtained from a regular run of FGBuster in pixel space, the $N_{side}$ = 512 maps were estimated by applying the W operator (eq. (3.21)) in harmonic space in order to have an optimal weighting of each frequency map by its associated beam.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Comparison of the pixel-based results using different FGBuster runs at different resolutions, and results using the CAB-SECRET maps. CAB-SECRET maps are used at $N_{side} = 64$ with a beam FWHM of 132 arcminutes. The fiducial angle ('Fid.') is zero for Phases 1, 2 and 3, and 0.3° for Phase 4. 'Est.' is the central value of the estimated $\hat{\beta}$ from the 100 simulations, the error bars are given by the standard deviation from the 100 simulations. The best results were achieved using FGBuster $N_{side} = 64$ for Phases 1 and 2, and FGBuster $N_{side} = 128$ for Phases 3 and 4. The $N_{side}$ 128 resolution FGBuster maps outperformed the highest resolution maps due to the additional noise complexity in the highest resolution maps.</figcaption></figure>

$N_{side} = 512$ maps perform less well for the stacking approach in particular. Other component separation techniques may have improvements at the $N_{side} = 512$ resolution in future work. The large beam of the CAB-SECRET maps reduces the ability of the peak-stacking method to estimate the final result with high accuracy. Histograms showing the results for the FGBuster Phase 1 and 2 at $N_{side} = 64$, and Phase 3 and 4 at $N_{side} = 128$ are presented in figure 10.

As shown in table 4, most of the constraining power from the stacking method comes from the extrema in the *E* maps, with a small boost from the extrema in the *T* maps. This is expected because the *E* peaks correspond to the *EB* signal from the power spectrum, whereas the *T* peaks correspond to the *TB* signal from the power spectrum. The justification for the higher signal-to-noise of *EB* with respect to *TB* can be derived from the Fisher information matrix [35, 99]. The consistency of $\beta$ estimates can be assessed by stacking only hot or only cold peaks and the comparison of the various data splits can help ascertain to what degree foregrounds may be affecting results, since foregrounds might be expected to change cold or hot peaks differently. Furthermore, the local nature of the estimator can be used to test the isotropy of the signal by targeting peaks in different regions of the sky.

## 4.2 Pipelines with control of instrumental polarisation angles

### 4.2.1 Template-based Minami-Komatsu estimator

As mentioned in section 3.2.1, the accuracy of the template-based MK estimator depends on our ability to model the *EB* correlation of Galactic foreground emission. Therefore, we focus our discussion on studying the impact that foreground modeling has on our measurement of CB. Here, we limit our analysis to s0, s1, d0, and d1 models and leave the study of further synchrotron and dust models to future work.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Results for 100 *LiteBIRD* simulation maps from the pixel-based analysis, generated using FGBuster. Details regarding each phase may be found in the text and numerical averages in table 4.</figcaption></figure>

<table><caption>**Table 4.** Comparison of the pixel-based results for the best data set from each phase (i.e., FGBuster, $N_{side}$ 64 for Phases 1 and 2, and FGBuster $N_{side}$ 128 for Phases 3 and 4), showing the average results from hot and cold spots for the *E* and *T* peaks (see eq. (4.3) for calculation). Results between each of the data splits should be consistent with one another. Inconsistencies would be an indication of effects due to uncontrolled systematics, such as foreground contamination.</caption><thead><tr><th rowspan="2">Phase</th><th><i>T</i> cold peaks</th><th><i>T</i> hot peaks</th><th><i>E</i> cold peaks</th><th><i>E</i> hot peaks</th><th>All</th></tr><tr><th><i>β</i> [deg]</th><th><i>β</i> [deg]</th><th><i>β</i> [deg]</th><th><i>β</i> [deg]</th><th><i>β</i> [deg]</th></tr></thead><tbody><tr><td>1</td><td>-0.024 ± 0.104</td><td>-0.008 ± 0.096</td><td>0.001 ± 0.046</td><td>-0.003 ± 0.048</td><td>-0.003 ± 0.030</td></tr><tr><td>2</td><td>-0.027 ± 0.121</td><td>-0.005 ± 0.110</td><td>-0.002 ± 0.057</td><td>-0.003 ± 0.060</td><td>-0.005 ± 0.037</td></tr><tr><td>3</td><td>0.000 ± 0.094</td><td>0.001 ± 0.103</td><td>0.011 ± 0.076</td><td>-0.005 ± 0.066</td><td>0.002 ± 0.040</td></tr><tr><td>4</td><td>0.285 ± 0.101</td><td>0.294 ± 0.101</td><td>0.285 ± 0.076</td><td>0.272 ± 0.063</td><td>0.281 ± 0.040</td></tr></tbody></table>

We start our discussion with Phase 2 simulations (where $\beta, \alpha_f = 0$) because this is the case that allows us to isolate the contribution of the foreground *EB*. Table 5 shows the mean CB and template amplitudes recovered with the MK estimator when the foreground *EB* correlation is ignored or modelled through different templates. The corresponding $\alpha_f$ are shown in figure 11. Error bars correspond to the 68% C.L. uncertainties calculated from the values' dispersion. As described in ref. [50], we expect foreground emission to bias our estimates to $\hat{\beta} \approx \beta^{\text{fid}} - \gamma_{\ell}$ and $\hat{\alpha}_f \approx \alpha_f^{\text{fid}} + \gamma_{\ell}$, where $\gamma_{\ell} \approx C_{\ell}^{EB,fgs}/(C_{\ell}^{EE,fgs} - C_{\ell}^{BB,fgs})$.
<table><caption>**Table 5.** Mean CB angle and template amplitudes recovered from Phase 2 simulations. Uncertainties are calculated as the standard deviation over 100 simulations (68% C.L.). At every row, we only fit for the amplitudes indicated in the first column. Note the sign change compared to figure 11, which shows $\delta_{\beta} = \langle 0 - \hat{\beta} \rangle = -\langle \hat{\beta} \rangle$. We use the Akaike information criterion (AIC) to compare the goodness of fit of the different models, evaluating $\Delta$AIC = AIC$_i$ - AIC$_{ref}$ with the $A^{sync} + A^{sync×dust} + A^{dust}$ (d1s1) case as reference.</caption><thead><tr><th>Foreground template</th><th>β [×10<sup>-2</sup>deg]</th><th>A<sup>sync</sup></th><th>A<sup>sync×dust</sup></th><th>A<sup>dust</sup></th><th>ΔAIC</th></tr></thead><tbody><tr><td>A<sup>x</sup> = 0</td><td>-1.12 ± 6.42</td><td>...</td><td>...</td><td>...</td><td>618</td></tr><tr><td>A<sup>dust</sup> (d0)</td><td>-0.29 ± 1.22</td><td>...</td><td>...</td><td>0.9997 ± 0.0006</td><td>343</td></tr><tr><td>A<sup>dust</sup> (d1)</td><td>-0.24 ± 1.61</td><td>...</td><td>...</td><td>1.0002 ± 0.0005</td><td>252</td></tr><tr><td>A<sup>sync</sup> + A<sup>dust</sup> (d0s0)</td><td>-0.35 ± 1.21</td><td>0.93 ± 0.05</td><td>...</td><td>0.9998 ± 0.0006</td><td>315</td></tr><tr><td>A<sup>sync</sup> + A<sup>dust</sup> (d1s1)</td><td>-0.16 ± 1.59</td><td>0.97 ± 0.06</td><td>...</td><td>1.0001 ± 0.0005</td><td>214</td></tr><tr><td>A<sup>sync</sup> + A<sup>sync×dust</sup> + A<sup>dust</sup> (d0s0)</td><td>-0.39 ± 1.22</td><td>0.94 ± 0.06</td><td>1.01 ± 0.02</td><td>0.9998 ± 0.0006</td><td>211</td></tr><tr><td>A<sup>sync</sup> + A<sup>sync×dust</sup> + A<sup>dust</sup> (d1s1)</td><td>-0.15 ± 1.66</td><td>0.99 ± 0.08</td><td>1.00 ± 0.02</td><td>1.0001 ± 0.0005</td><td>0</td></tr></tbody></table>

Note that $\hat{\alpha}_f + \hat{\beta}$ remains unaffected, since we are trying to use foregrounds to break the degeneracy between both rotations. Hence, the small $C_{\ell}^{EB,fgs} > 0$ found in the s0, s1, d0, and d1 models produces the $\hat{\beta} < \beta^{fid}$ and $\hat{\alpha} > \alpha^{fid}$ biases seen in figure 11. Similar biases were also found in refs. [50, 55].

When neglected ($A^{dust}$, $A^{sync}$, $A^{sync×dust} = 0$ in eq. (3.10) and "$A_x = 0$" entries in table 5 and figure 11), the EB correlation in d1 produces a bias of $\delta_{\beta} = \langle \beta^{fid} - \hat{\beta} \rangle = 0.011^{\circ}$ ($0.17\sigma$ bias). This bias is lower than the $\delta_{\beta} \approx 0.05^{\circ}$ to $0.15^{\circ}$ biases reported in almost-full-sky analyses of Planck data<fnref n="12" /> [50], partially because of the simplified setup of our simulation framework, but also thanks to the configuration of LiteBIRD itself. The numerous cross-correlations between dust- and synchrotron-dominated bands possible within LiteBIRD's 22 frequency channels help suppress the impact of dust EB [51, 52]. Neglecting dust EB also biases $\alpha_f$ estimates [55, 91]. Due to their higher noise levels and lower resolutions [14], $\alpha_f$ estimation on the LFT's lower frequencies relies mostly on foreground-contaminated scales, leading to the larger biases and uncertainties seen in those bands [91].

Providing a dust template ($A^{dust} \neq 0$ while $A^{sync}$, $A^{sync×dust} = 0$ in eq. (3.10) and "$A^{dust}$ (d1)” entries in table 5 and figure 11) is enough to reduce the dust-EB bias to $\delta_{\beta} = 0.0024^{\circ}$ and bring both the $\beta$ and $\alpha_f$ angles of dust-dominated frequency bands in to agreement with zero at the 0.12 to 0.26$\sigma$ level. As discussed in ref. [55], the mode-by-mode subtraction of the foreground signal produced by the template also leads to an overall reduction of uncertainties. For the lower frequencies of the LFT, $\alpha_f$ estimation remains dominated by the instrumental configuration, since the dust template does not provide additional constraining power at the frequencies where dust is subdominant.

We expect synchrotron radiation to have a negligible impact on the MK estimator [51, 52] since no significant synchrotron EB has been found to date [86, 87]. This is confirmed by

<footnote n="12">¹²For almost-full-sky analyses ($f_{sky} \gtrsim 90\%$), we expect the local misalignments between different dust filaments and the Galactic magnetic field to average out along the line of sight and give a small global EB correlation [53]. Depending on the dust model, a $\delta_{\beta} \approx 0.05^{\circ}$ to $0.15^{\circ}$ bias is expected for this $f_{sky}$ [50, 51]. Far from the Galactic plane, a larger spatial and angular coherence is seen in the misalignment of dust filaments and magnetic fields [54], leading to $\delta_{\beta} \approx 0.55^{\circ}$ biases for $f_{sky} \approx 60\%$.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Average difference between the input and estimated angles from Phase 2 simulations when the foreground EB correlation is ignored or modelled with different templates. Uncertainties show the $1\sigma$ dispersion over simulations. Providing a dust template suffices to obtain an unbiased estimate of $\beta$, and MFT's and HFT's $\alpha$ at the $0.12\sigma$ to $0.26\sigma$ level (depending on the frequency band). Extra information from a synchrotron template is needed to correctly recover $\alpha$ at the lower LFT frequencies. The EB correlation in d0 is a good enough approximation of that in d1 to effectively correct the dust EB bias. In contrast, s0 proves an insufficient description of s1, limited in its improvement on the calibration of $\alpha$.</figcaption></figure>

the small impact including a synchrotron template ($\mathcal{A}^{\text{sync}}$, $\mathcal{A}^{\text{dust}} \neq 0$ while $\mathcal{A}^{\text{sync}\times\text{dust}} = 0$ in eq. (3.10) and "$\mathcal{A}^{\text{sync}} + \mathcal{A}^{\text{dust}} (\text{d1s1})$" entries in table 5 and figure 11) has on the CB measurement once we include a dust template. Nevertheless, the additional information provided by the synchrotron template helps to characterise the large angular scales of the observed signal, allowing a better estimation of $\alpha_{40}$, $\alpha_{50}$, $\alpha_{60}$, $\alpha_{68\text{a}}$, and $\alpha_{68\text{b}}$ despite the LFT's limited angular resolution. Likewise, providing additional information on the correlation between dust and synchrotron emission ($\mathcal{A}^{\text{sync}}$, $\mathcal{A}^{\text{dust}}$, $\mathcal{A}^{\text{sync}\times\text{dust}} \neq 0$ in eq. (3.10) and "$\mathcal{A}^{\text{sync}} + \mathcal{A}^{\text{sync}\times\text{dust}} + \mathcal{A}^{\text{dust}} (\text{d1s1})$" entries in table 5 and figure 11) further improves the estimate of $\alpha_f$ in the LFT's low-frequency bands. With the observed signal and $\alpha_f$ better constrained, more information about $\beta$ can be gathered from low-frequency bands. Thus, providing a synchrotron template improves the overall performance of the estimator and reduces the foreground-induced bias to $\delta_{\beta} = 0.0015^\circ$.

An accurate estimation of $\beta \approx 0.3^\circ$ is still possible under a small mismodelling of foreground emission, like assuming d0s0 templates to describe an d1s1 sky. Models d0 and d1 are based on the same dust template derived from Planck 2015 data, but scaled with different SEDs: a fixed $\beta_{\text{d}} = 1.54$ and $T_{\text{d}} = 20$ K for d0; and the spatially varying $T_{\text{d}}$ and
$\beta_d$ obtained by Commander from the analysis of Planck 2015 data [100] for d1. Likewise, s0 and s1 share the synchrotron template derived from Haslam and WMAP data [101, 102], but are scaled with different SEDs: a fixed $\beta_s = -3.0$ for s0; and the spatially varying $\beta_s$ obtained by [103] from the analysis of Haslam and WMAP data for s1. See refs. [62, 90] for further detail. At the angular scales used in our analysis, only 5% to 15% differences are found between the EB spectra of d0 and d1, and s0 and s1.

Qualitatively, d0 offers a good enough description of d1 to satisfactorily correct its EB correlation as the minor increase in the $\delta_\beta$ bias and the good agreement of $A^{\text{dust}}$ with unity show in table 5 ("$A^{\text{dust}}$ (d0)" entries). On the contrary, in our joint analysis of dust and synchrotron templates, s0 proves insufficient to describe s1 ("$A^{\text{sync}} + A^{\text{dust}}$ (d0s0)" and "$A^{\text{sync}} + A^{\text{sync}\times\text{dust}} + A^{\text{dust}}$ (d0s0)" entries in table 5 and figure 11): modeling s1 with s0 significantly increases the $\delta_\beta$ bias, gives $A^{\text{sync}}$ values incompatible with unity within $1\sigma$, and no longer improves the estimation of $\alpha_f$ at synchrotron-dominated frequencies. In addition to increasing the bias on $\beta$ and $\alpha_f$, assuming simplified synchrotron and dust models in the calculation of the covariance matrix also leads to roughly a 25% underestimation of uncertainties.

Since it is based on a maximum-likelihood algorithm, our MK estimator allows the assessment of the goodness of fit and the comparison of different foreground models through statistics like the Akaike information criterion (AIC). We complement our discussion by showing the $\Delta$AIC values obtained with different foreground templates in table 5. The recovered $\Delta$AIC reflect the inadequacy of neglecting foreground EB and confirm $A^{\text{sync}} + A^{\text{sync}\times\text{dust}} + A^{\text{dust}}$ as the best description of the foreground emission present in the simulations.

As figure 12 and table 6 show, these conclusions qualitatively hold when moving to the more general scenario of Phase 4 simulations, where $\beta$ and $\alpha_f$ take non-zero values. In this case, the observed signal is more complex than in Phase 2 simulations, since the E and B modes of Galactic foregrounds are rotated into each other. The complexity of the covariance matrix describing such a signal also increases, gaining numerous terms that depend on the cross-correlations between observations and foreground templates (see appendix A). In this more realistic scenario, providing an accurate foreground template is crucial. Without one, the estimated $\beta$ would suffer a bias of $\delta_\beta = 0.1164^\circ$ ($1.8\sigma$ bias).

With a foreground template, the MK estimator accurately measures $\beta \approx 0.3^\circ$ with a small bias of $\delta_\beta \approx 0.008^\circ$ (40% of the statistical uncertainty). Additionally, we recover a stable $\delta_\beta$ value for all the sky models considered, since now the inaccuracy of the approximations made in the likelihood dominates over the effects of foreground mismodelling. Those approximations consist of assuming that the statistics of angular power spectra follow a Gaussian distribution and using an approximate covariance matrix derived from the observed signal on a realisation-by-realisation basis. We further simplify the covariance matrix by neglecting $\ell$-to-$\ell'$ correlations; assuming that the spherical harmonic coefficients of CMB, noise, and foreground signals are Gaussian and isotropic; and discarding some of the more noisy and unstable EB terms of the observed signal. While these approximations proved sufficient for the simpler framework of Phase 2, they fall short for Phase 4 simulations. For instance, the $d^T C^{-1} d$ product of $C_\ell$ data vector and covariance matrix that followed the expected $\chi^2$ distribution for Phase 2 simulations, no longer does so for Phase 4. In practice, this means that, although sufficient to provide an unbiased estimate of $\beta$, the implementation of the MK
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Same as figure 11, but for Phase 4 simulations. For Phase 4, providing a dust template allows an unbiased estimate of $\beta$, and MFT's and HFT's $\alpha$ at the $0.19\sigma$ to $0.29\sigma$ level. Similar conclusions on the need for a synchrotron template to accurately calibrate LFT's $\alpha$ and the (in)adequacy of (s0) d0 models are obtained for Phase 4.</figcaption></figure>

<table><caption>**Table 6.** Same as table 5 but for Phase 4 simulations. Although sufficient to provide an unbiased CB estimate, the current MK estimator is not suitable for model comparison through statistics like the Akaike information criterion in the more complex scenario of Phase 4.</caption><thead><tr><th>Foreground template</th><th>β [×10<sup>-2</sup> deg]</th><th>A<sup>sync</sup></th><th>A<sup>sync×dust</sup></th><th>A<sup>dust</sup></th></tr></thead><tbody><tr><td>A<sup>x</sup> = 0</td><td>18.36 ± 6.62</td><td>...</td><td>...</td><td>...</td></tr><tr><td>A<sup>dust</sup> (d0)</td><td>29.20 ± 1.81</td><td>...</td><td>...</td><td>1.002 ± 0.003</td></tr><tr><td>A<sup>dust</sup> (d1)</td><td>29.25 ± 1.98</td><td>...</td><td>...</td><td>0.998 ± 0.002</td></tr><tr><td>A<sup>sync</sup> + A<sup>dust</sup> (d0s0)</td><td>29.10 ± 1.88</td><td>0.92 ± 0.06</td><td>...</td><td>0.999 ± 0.004</td></tr><tr><td>A<sup>sync</sup> + A<sup>dust</sup> (d1s1)</td><td>29.28 ± 2.11</td><td>0.96 ± 0.09</td><td>...</td><td>0.999 ± 0.001</td></tr><tr><td>A<sup>sync</sup> + A<sup>sync×dust</sup> + A<sup>dust</sup> (d0s0)</td><td>29.23 ± 1.46</td><td>0.92 ± 0.07</td><td>0.99 ± 0.03</td><td>0.999 ± 0.004</td></tr><tr><td>A<sup>sync</sup> + A<sup>sync×dust</sup> + A<sup>dust</sup> (d1s1)</td><td>29.27 ± 2.24</td><td>0.96 ± 0.13</td><td>0.99 ± 0.04</td><td>0.999 ± 0.002</td></tr></tbody></table>

estimator presented here is not suitable for model comparison through statistical estimators like AIC in the high signal-to-noise regime offered by *LiteBIRD*. Therefore, we must explore extensions of the estimator beyond the Gaussian likelihood (see, e.g., ref. [104]) and a more robust calculation of its covariance matrix if we want to reliably assess the goodness of fit of different foreground models once *LiteBIRD* data arrive. We leave that study to future work.
<table><caption>**Table 7.** Biases and uncertainties for each phase obtained using either CAB-SeCRET or βCAB-SeCRET. The biases are expressed in terms of the uncertainty on β. σ<sub>β</sub> and σ<sup>MC</sup><sub>β</sub> represent, respectively, the standard deviation calculated from simulations and from the MCMC chains obtained by minimizing eq. (3.15).</caption><thead><tr><th>Phase</th><th>Method</th><th>Prior</th><th>N<sup>B</sup><sub>side</sub></th><th>δ<sub>β</sub>/σ<sub>β</sub></th><th>σ<sub>β</sub> [deg]</th><th>δ<sup>MC</sup><sub>β</sub>/σ<sup>MC</sup><sub>β</sub></th><th>σ<sup>MC</sup><sub>β</sub> [deg]</th></tr></thead><tbody><tr><td>1</td><td>CAB-SeCRET</td><td>Π(α<sub>f</sub>, β = 0)</td><td>0</td><td>-0.104</td><td>0.0336</td><td>-0.124</td><td>0.0282</td></tr><tr><td>2</td><td>CAB-SeCRET</td><td>Π(α<sub>f</sub>, β = 0)</td><td>0</td><td>-2.378</td><td>0.0357</td><td>-2.891</td><td>0.0294</td></tr><tr><td>2</td><td>CAB-SeCRET</td><td>Π(α<sub>f</sub>, β = 0)</td><td>8</td><td>0.420</td><td>0.0341</td><td>0.407</td><td>0.0352</td></tr><tr><td>2</td><td>CAB-SeCRET</td><td>Π(α<sub>f</sub>, β = 0)</td><td>64</td><td>-0.218</td><td>0.0840</td><td>-0.279</td><td>0.0656</td></tr><tr><td>3</td><td>CAB-SeCRET</td><td>Π(α<sub>f</sub>, β = 0)</td><td>8</td><td>0.483</td><td>0.0349</td><td>0.479</td><td>0.0352</td></tr><tr><td>4</td><td>CAB-SeCRET</td><td>Π(α<sub>f</sub>, β)</td><td>8</td><td>8.414</td><td>0.0373</td><td>10.056</td><td>0.0312</td></tr><tr><td>4</td><td>βCAB-SeCRET</td><td>Π(α<sub>f</sub>, β)</td><td>8</td><td>-0.160</td><td>0.0359</td><td>-0.167</td><td>0.0345</td></tr></tbody></table>

Finally, these results demonstrate that *LiteBIRD* is capable of self-calibrating instrumental α<sub>f</sub> and measuring isotropic CB at the same time, as long as a reasonable foreground model is provided. With the help of synchrotron and dust templates, we can recover unbiased estimates of α<sub>f</sub> with uncertainties of around 2 arcmin for all frequency bands, fulfilling the requirements needed to achieve an unbiased measurement of the tensor-to-scalar ratio with *LiteBIRD* [60]. This study extends previous results presented in refs. [46, 48, 64, 91].

### 4.2.2 CAB-SeCRET

In this section, we discuss the results obtained with CAB-SeCRET and βCAB-SeCRET, focusing on the subtleties intrinsic to this pipeline. First, we present the results obtained for each phase using CAB-SeCRET and discuss the impact of mismodelled spectral parameters and biased priors on α<sub>f</sub>. Then, we compare the Phase 4 results recovered with CAB-SeCRET and βCAB-SeCRET.

**CAB-SeCRET results.** This pipeline takes into account the effect of nonzero α<sub>f</sub> in the component-separation analysis, as inaccuracies in α<sub>f</sub> calibration can introduce substantial biases in CB detection. Figure 13 shows the average difference between the input and recovered β and α<sub>f</sub> obtained with CAB-SeCRET for each phase. We apply the Π(α<sub>f</sub>, β = 0) priors for Phases 1 to 3 and switch to Π(α<sub>f</sub>, β) for Phase 4. Had we employed Π(α<sub>f</sub>, β) for Phases 1 to 3 as well, we would have observed a slight increase in σ<sub>β</sub>. However, this increase is deemed insignificant upon comparing the uncertainty observed in Phase 4 to that of the preceding Phases, as illustrated in table 7.

Subsequently, we examine the outcomes for each phase, in order of ascending complexity.

* **Phase 1:** for Phase 1, we assume that the *B* parameters are uniform across the entire observable sky after applying the *f*<sub>sky</sub> = 60% Planck Galactic mask. This modelling approach reproduces that of the simulated foregrounds, ensuring that we do not obtain biased results due to incorrect modeling of the sky. Our findings indicate that the α<sub>f</sub> overall align with the values injected in the simulation, i.e., α<sub>f</sub> = 0, as shown in
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Average difference between the input and estimated angles obtained with CAB-SECRET. Note that $\alpha_f$ are calculated at the component-separation step, while $\beta$ is calculated from the cleaned CMB map. Blue diamonds, orange pluses, and green circles correspond, respectively, to Phase 1, Phase 2, and Phase 3 results obtained using the $\Pi(\alpha_f, \beta = 0)$ priors. Red squares show the Phase 4 results when we apply the $\Pi(\alpha_f, \beta)$ priors. Uncertainties correspond to the standard deviation of the simulations (68% C.L.). Note the difference in the y-axis scale between the top and bottom panels.</figcaption></figure>

figure 13. Furthermore, we demonstrate that the $\beta$ left in the clean CMB map is consistent with zero, with an uncertainty of $\sigma_{\beta} \simeq 0.034^{\circ}$ and a bias of $\delta_{\beta} \approx -0.1\sigma_{\beta}$.

* **Phase 2:** for this simulation set, the assumption of uniform $\mathcal{B}$ parameters no longer holds. Relaxing this condition implies the addition of more free parameters, increasing the uncertainty of all model parameters, including $\beta$. Thus, we need to find a trade-off between a tolerable $\delta_{\beta}$ and the $\sigma_{\beta}$ necessary to achieve a detection if $\beta$ is of the order of $0.3^{\circ}$. We consider three different cases.

i. As for Phase 1, we assume that $\mathcal{B}$ parameters are uniform, meaning the resolution of $\mathcal{B}$ is $N_{side}^{\mathcal{B}} = 0$.

ii. We allow $\mathcal{B}$ to vary from pixel to pixel ($N_{side}^{\mathcal{B}} = 64$). This case aligns the most with the fiducial foreground model, hence, we expect negligible biases.

iii. An intermediate case, where $\mathcal{B}$ are assumed to be uniform in superpixels of $N_{side}^{\mathcal{B}} = 8$. Here, we incorporate spatial variations of $\mathcal{B}$ but with fewer parameters, approximately a factor of $3 \times 64$ fewer than in the previous one.

As anticipated, the value of $\sigma_{\beta}^{\text{MC}}$ increases with the number of $\mathcal{B}$ parameters ($N_{side}^{\mathcal{B}}$). Explaining the value of $\sigma_{\beta}$ is more challenging, since it accounts for both the dispersion due to a larger number of $\mathcal{B}$ parameters and the bias resulting from imperfect modelling
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Left: distribution difference between the estimated $\alpha_{\text{out}}$ for all *LiteBIRD* channels and the input $\alpha_f$ (blue line histogram), and the input $\alpha_f$ plus the injected $\beta$ angle (orange filled histogram) when we fit using the $\Pi(\alpha_f, \beta = 0)$ priors. Right: distribution of the estimated $\beta$ when we apply the priors assuming $\beta = 0$ (blue line histogram) and with no constraints on $\beta$ (orange filled histogram). The dashed grey line indicates the input $\beta$ value.</figcaption></figure>

of the foreground emission. Regarding the bias, if $N_{\text{side}}^{\mathcal{B}} = 0$ like in the previous Phase, we recover a biased value of $\delta_{\beta} = -2.379\sigma_{\beta}$ due to imperfect modelling of contaminants. When including all possible parameters at the current resolution, i.e., $N_{\text{side}}^{\mathcal{B}} = 64$, allowing the model to capture the complexity of the spatial variations of spectral indices, the uncertainty increases to $\sigma_{\beta} = 0.084^{\circ}$, corresponding to a $\beta$ detection of $3.6\sigma$, similar to current constraints. With $N_{\text{side}}^{\mathcal{B}} = 8$, we find a compromise between bias, $\delta_{\beta} = 0.420\sigma_{\beta}$, and competitive uncertainty, $\sigma_{\beta} = 0.0341^{\circ}$. We use $N_{\text{side}}^{\mathcal{B}} = 8$ in the following phases.

* **Phase 3:** figure 13 demonstrates that the estimates for both $\alpha_f$ and $\beta$ match the input values when non-zero $\alpha_f$ are introduced in the simulations. Furthermore, we observe that $\delta_{\beta}$ and $\sigma_{\beta}$ from Phase 3 closely resemble the values from Phase 2 ($N_{\text{side}}^{\mathcal{B}} = 8$), confirming a robust marginalisation over $\alpha_f$.

* **Phase 4:** in this set of simulations, $\beta$ is no longer zero, rendering the method to estimate $\alpha_f$ from ref. [91] inaccurate as it estimates the combined angle $\tilde{\alpha}_f = \alpha_f + \beta$. Since $\beta$ represents a common rotation across all frequency channels, $\tilde{\alpha}_f$ also provides a reasonable fit to the data. Consequently, the recovered $\alpha_f$ are biased. This is evident from the left panel of figure 14, where we compare the distribution of the difference between the recovered $\tilde{\alpha}_f$ and, firstly, the injected $\alpha_f^{\text{in}}$ in the simulations, and secondly, $\alpha_f^{\text{in}} + \beta$. The former distribution peaks around $0.3^{\circ}$, corresponding to the input $\beta$ value, while the latter peaks around zero.

To address this issue, one can use prior information regarding $\alpha_f$ calibration. In particular, we use the $\mathcal{A}^{\text{sync}} + \mathcal{A}^{\text{sync}\times\text{dust}} + \mathcal{A}^{\text{dust}}$ (d1s1) case from section 4.2.1 in the following. It is important to highlight that we do not explicitly incorporate any $\beta$ constraints from the MK estimator. We only make use of the information on $\alpha_f$ calibration that could otherwise be provided by any other absolute calibration method.
<table><caption>**Table 8.** Biases and uncertainties from Phase 4 obtained using the template-based MK ($A^{sync} + A^{sync \times dust} + A^{dust}$ with d1s1), CAB-SeCRET, and βCAB-SeCRET methods. The biases are expressed in terms of the uncertainty on $\beta$. The quantity $\sigma_{\beta}$ represents the standard deviation obtained from simulations.</caption><thead><tr><th>Method</th><th>δ<sub>β</sub>/σ<sub>β</sub></th><th>σ<sub>β</sub> [deg]</th></tr></thead><tbody><tr><td>Template-based MK</td><td>-0.262</td><td>0.0214</td></tr><tr><td>CAB-SeCRET</td><td>0.160</td><td>0.0359</td></tr><tr><td>βCAB-SeCRET</td><td>-0.037</td><td>0.0617</td></tr></tbody></table>

Additionally, in this approach, we simultaneously marginalise over the foreground model parameters and $\alpha_f$. Thus, we alleviate any potential biases associated with the mismodelling of the foreground templates in the MK estimator.

With these priors, we obtain $\alpha_f$ and $\beta$ values that are consistent with the input values, as depicted in figure 13, and right panel of figure 14. Although the uncertainty in the recovered $\alpha_f$ increases due to the less stringent priors, we find that the uncertainty in the recovered $\beta$ is similar to that of the other phases. This is because $\beta$ is derived from the cleaned CMB, whose noise is minimally affected by changes in the $\alpha_f$ prior. However, we observe a slight reduction in the $\delta_{\beta}$ bias, despite the astrophysical model being incorrect, as in the other Phases 2 and 3.

**Comparison of CAB-SeCRET and βCAB-SeCRET.** In this section, we compare the results from Phase 4 obtained by applying the CAB-SeCRET and βCAB-SeCRET. Notice that, in both cases, we use the results from the template-based MK $A^{sync} + A^{sync \times dust} + A^{dust}$ (d1s1) estimator as prior information ($\Pi(\alpha_f, \beta)$). In CAB-SeCRET, we use only the information from $\alpha_f$, while in βCAB-SeCRET we also incorporate the information from $\beta$. The results are summarised in table 8 and illustrated in figure 15, where we also include the prior estimates for comparison. When comparing the results from CAB-SeCRET with respect to βCAB-SeCRET, we observe an increase in the uncertainty on $\beta$ and a reduction in the $\delta_{\beta}$ bias.

The primary distinction between CAB-SeCRET and βCAB-SeCRET lies in their treatment of $\beta$: in the former, $\beta$ is conditioned on both $\alpha_f$ and the astrophysical component parameters, while in the latter, $\beta$ is marginalised over them. Consequently, the results from βCAB-SeCRET are more robust and unbiased, as demonstrated in table 8. Furthermore, it is noteworthy that although the uncertainty in $\beta$ obtained with βCAB-SeCRET is larger, it remains competitive, especially if $\beta \approx 0.3^{\circ}$, as it would result in a detection of $4.9\sigma$.

Similar conclusions arise from comparing the results of CAB-SeCRET and βCAB-SeCRET with those from the template-based MK estimator used as prior information. Both CAB-SeCRET and βCAB-SeCRET simultaneously marginalise over astrophysical model parameters and $\alpha_f$, unlike the template-based MK estimator. This marginalisation results in smaller biases but larger uncertainties. Consequently, the two pipelines can complement each other. For instance, one can utilise the results from the template-based MK estimator to constrain the $C$ parameters and recover $\beta$ with βCAB-SeCRET or CAB-SeCRET, and then employ the foreground templates derived from the latter in the template-based MK estimator to mitigate biases arising from incorrect foreground templates.
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Left: distribution of the estimated $\beta$ from Phase 4 simulations using three methods: CAB-SECRET; $\beta$CAB-SECRET; and the template-based MK results ($A^{\text{sync}} + A^{\text{sync}\times\text{dust}} + A^{\text{dust}}$ with d1s1) used as prior information in $\beta$CAB-SECRET. Right: average difference between input and estimated $\beta$.</figcaption></figure>

### 4.2.3 Pipeline J23

In this section, we discuss the results obtained by the J23 pipeline, going phase by phase and discussing the possible shortcomings and extensions to the pipeline in view of its response to the different levels of complexity of each phase. All error bars discussed here correspond to the standard deviation over the set of simulations: $\sigma_x$ (68% C.L.) for parameter $x$.

* **Phase 1:** the foreground cleaning in this Phase should be close to optimal, as the simulations correspond to the data model used in the component separation (eq. (3.19)). Furthermore, the absence of $\alpha_f$ miscalibration should lead to a reduced $\sigma_\beta$ compared to phases with miscalibration.

Figure 16 shows the results of the MCMC sampling of the generalised spectral likelihood (eq. (3.20)) for one given simulation, corresponding to step 1 in section 3.2.3. The orange contours correspond to the 68% and 95% contours from the MCMC sampling, while the purple dashed ones represent the Gaussian priors from calibration. Note that, as explained in section 3.2.3, the priors’ central values are chosen at random for each simulation, hence they are not centred at the true values of parameters denoted by the dotted grey lines. We only display $\alpha_f$ of three frequency channels, $\alpha_{\text{LFT,68a}}$, $\alpha_{\text{MFT,100}}$, and $\alpha_{\text{HFT,235}}$, as showing the full triangle plot with 22 angles would lack clarity. We can see that the foreground parameters are estimated correctly and seem uncorrelated to $\alpha_f$, as expected from ref. [95]. Thanks to the combination of prior information between frequencies and the internal calibration against the foregrounds SED, the statistical error bars associated with $\alpha_f$ estimates are smaller than the initial precision of the priors (which are $\sigma_{\alpha_{\text{LFT,68a}}} = 0.7^\circ$, $\sigma_{\alpha_{\text{MFT,100}}} = 0.13^\circ$, and $\sigma_{\alpha_{\text{HFT,235}}} = 0.21^\circ$) and best-fit values are closer to the true input angle. Although off-centred priors can bias $\alpha_f$, the mean bias on the absolute offset will tend to average out to zero thanks to the randomness of the priors’ centres. This can already be partially seen in the figure, with the first angle being biased positively while the two others are negatively biased.
Pipeline J23

<figure><img src="image_19.png" /></figure>

**Figure 16.** MCMC sampling of the generalised spectral likelihood (eq. (3.20)) for one Phase 1 simulation for pipeline J23. Out of the total 22 $\alpha_f$, only $\alpha_{LFT,68a}$, $\alpha_{MFT,100}$, and $\alpha_{HFT,235}$ are displayed for clarity. Orange contours represent the 68% and 95% contours. Purple dashed contours represent the Gaussian priors applied to the angles, as explained in section 3.2.3 (Step 1) the central value of the priors are drawn randomly for each simulation. True input values are displayed as the grey dotted lines.

Furthermore, the case presented in figure 16 only represents one simulation and thus does not encompass the full posterior and, in particular, the variance associated with the random prior centres.

Therefore, we must look at the results over the 100 simulations. Figure 17 shows the results for the estimation of $\beta$ and $\alpha_f$. Central values correspond to the average over simulations of the estimated parameters, and error bars to $\sigma_x$. We can see that all $\alpha_f$ are retrieved without significant bias, with the error bars on $\alpha_f$ ranging from 0.028° to 0.085°. This leads to an estimation of $\beta = -0.002° \pm 0.042°$, hence demonstrating
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Average difference between the input and estimated angles over Phase 1 simulations obtained with pipeline J23. Here, $\alpha_f$ are estimated in the component-separation step (section 3.2.3, Step 1), while $\beta$ is obtained using the cosmological likelihood with the angular power spectra of the foreground-cleaned CMB map (section 3.2.3, Step 5). Error bars represent the standard deviation of the values retrieved over the 100 simulations (68% C.L.).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 18.** From left to right, distribution of Phase 1 results for $\beta_d$, $T_d$, and $\beta_s$ for pipeline J23. Grey dashed lines represent the true input parameters used in the simulations.</figcaption></figure>

that, in this set of simple simulations, the J23 pipeline can correctly retrieve CB and $\alpha_f$. The histograms of best-fit values shown in figure 18 demonstrate that foreground parameters are also retrieved correctly, with average values of $\beta_d = 1.5401 \pm 0.0028$, $T_d = 19.99 \pm 0.10$, and $\beta_s = -3.0001 \pm 0.0039$.

* **Phase 2:** to handle the more complex case of spatially varying spectral parameters, CAB-SECRET and $\beta$CAB-SECRET allow for parameters to vary across superpixels during component separation. Contrary to these other pipelines, the current implementation of J23 does not have yet a way to handle cases with spatially varying spectral parameters and is therefore more susceptible to the biases caused by more complex foregrounds. Ref. [95] explored complex foregrounds (including d1s1 and d7s3,$^{13}$), but the smaller

<footnote>$^{13}$d7 uses a non-parametric dust frequency scaling based on dust population modelling [105]. s3 uses a power law with a curved index for the frequency scaling, where the index map is the same as for s1.</footnote>
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Similar to figure 17, but comparing Phase 2 with the standard $f_{\text{sky}} = 60\%$ mask (in blue) to Phase 2 with a larger mask $f_{\text{sky}} = 40\%$ (in orange) for pipeline J23.</figcaption></figure>

survey footprint of a Simons Observatory SAT-like survey ($f_{\text{sky}} \sim 10\%$) limits the impact of spatially varying foregrounds. This work is therefore the first to test this method on the bigger sky fraction typical of space missions.

Figure 19 summarises the Phase 2 results, showing that $\alpha_f$ at low frequencies are biased. In particular, the estimation of $\alpha_{40}$, the lowest frequency band, is biased by almost $5\sigma$. The bias on $\alpha_f$ reduces as the frequency increases, potentially indicating that this method is sensitive to spatially varying synchrotron emission. However, the estimation of $\alpha_f$ without a significant bias at middle and high frequencies partially compensates for the low-frequency bias, leading to an estimation of $\beta = -0.033^\circ \pm 0.047^\circ$. Hence, even without a way to handle complex foregrounds, the bias is below $1\sigma$ at $\delta_\beta = 0.708 \sigma_\beta$.

To confirm that the spatial variability of spectral parameters is indeed the reason for this bias, we perform a second analysis on Phase 2 simulations. In this case, we use a larger mask corresponding to $f_{\text{sky}} = 40\%$, since the smaller sky fraction should limit the impact of the foregrounds' spatial variation. The results are displayed in orange in figure 19 and in table 9. The bias on low-frequency $\alpha_f$ is indeed greatly reduced, with the bias on $\alpha_{40}$ now being $\lesssim 1\sigma$, and resulting in a CB estimate of $\beta = -0.010^\circ \pm 0.055^\circ$. The smaller sky fraction leads to a slight increase in $\sigma_\beta$ of about 15% and a drop in the bias to $\delta_\beta = 0.179 \sigma_\beta$.

Reducing the sky fraction to deal with complex foregrounds gives us insight into our analysis, but it should be treated as a diagnosis tool rather than an analysis tool. While having statistically consistent results on $\beta$ between $f_{\text{sky}} = 60\%$ and $f_{\text{sky}} = 40\%$ is promising, the large variation observed in the estimated $\alpha_f$ is a clear indication that complex foregrounds need to be dealt with. Possible extensions to the framework include the multi-patch approach [14, 106], which allows for spatial variability of foreground parameters, or foreground marginalisation [106, 107], using the spectra of the foreground maps estimated by the component separation in the cosmological likelihood. These extensions are beyond the scope of this paper and are left for future work.

For completeness, we also performed the analysis at $f_{\text{sky}} = 40\%$ on Phase 1 simulations, including the results in table 9. The smaller sky fraction leads again to an
<table><caption>**Table 9.** Summary of the results on $\beta$ obtained by pipeline J23 across the different phases and sky fractions.</caption><thead><tr><th>Phase</th><th>f<sub>sky</sub></th><th>β [deg]</th><th>δ<sub>β</sub>/σ<sub>β</sub></th></tr></thead><tbody><tr><td>1</td><td>0.6</td><td>-0.002 ± 0.042</td><td>0.049</td></tr><tr><td>1</td><td>0.4</td><td>-0.003 ± 0.051</td><td>0.051</td></tr><tr><td>2</td><td>0.6</td><td>-0.033 ± 0.047</td><td>0.708</td></tr><tr><td>2</td><td>0.4</td><td>-0.010 ± 0.055</td><td>0.179</td></tr><tr><td>3</td><td>0.6</td><td>-0.041 ± 0.048</td><td>0.849</td></tr><tr><td>3</td><td>0.4</td><td>-0.018 ± 0.056</td><td>0.320</td></tr><tr><td>4</td><td>0.6</td><td>0.258 ± 0.049</td><td>0.852</td></tr><tr><td>4</td><td>0.4</td><td>0.282 ± 0.056</td><td>0.322</td></tr></tbody></table>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** Similar to figure 19, comparing results of Phase 3 with $f_{\text{sky}} = 60\%$ (blue) and $f_{\text{sky}} = 40\%$ (orange) for pipeline J23.</figcaption></figure>

increase of around 20% in $\sigma_{\beta}$, but there is no significant difference in the bias thanks to the absence of spatial variability in these simulations.

* **Phase 3:** for this phase we proceed similarly, with a version of the analysis performed with the fiducial mask at $f_{\text{sky}} = 60\%$ and another with $f_{\text{sky}} = 40\%$. Results are displayed in figure 20 in blue and orange, respectively, and are summarized in table 9. Note that, since figure 20 shows the mean and standard deviation of the difference between true and estimated angles. Overall, this results in an increase of 3.3% of the error bar on $\beta$ compared to the Phase 2 result for the $f_{\text{sky}} = 60\%$ case, and 1.7% for the $f_{\text{sky}} = 40\%$ cases. The retrieved mean value of $\beta$ is also affected, increasing by 22% ($f_{\text{sky}} = 60\%$) and 77% ($f_{\text{sky}} = 40\%$) going from Phase 2 to 3. Results are still within the $1\sigma$ limit, with relative bias going to $\delta_{\beta} = 0.849\sigma_{\beta}$ ($f_{\text{sky}} = 60\%$) and $\delta_{\beta} = 0.320\sigma_{\beta}$ ($f_{\text{sky}} = 40\%$). Hence, in the presence of additional systematic effects, pipeline J23 seems to give robust results. However, the discussion in Phase 2 about requiring the extension of the pipeline to handle complex foregrounds holds here.
<figure><img src="image_24.png" /><figcaption>**Figure 21.** Similar to figures 19 and 20, comparing results of Phase 4 with $f_{\text{sky}} = 60\%$ (blue) and $f_{\text{sky}} = 40\%$ (orange) for pipeline J23.</figcaption></figure>

* **Phase 4:** finally, for Phase 4, we display the results in figure 21 and in table 9. Results are similar to Phase 3, with no significant impact on the estimation of $\alpha_f$. CB is retrieved with $\beta = 0.258^\circ \pm 0.049^\circ$ ($f_{\text{sky}} = 60\%$) and $\beta = 0.282^\circ \pm 0.056^\circ$ ($f_{\text{sky}} = 40\%$), corresponding to $\delta_\beta = 0.852 \sigma_\beta$ and $\delta_\beta = 0.322 \sigma_\beta$, respectively, and demonstrating the capacity of pipeline J23 to detect $\beta = 0.3^\circ$ at the $5\sigma$ level.

These results show that *LiteBIRD*, with the use of pipeline J23, is capable of measuring $\beta$ with a precision of about $0.05^\circ$. This is possible thanks to our knowledge of the foreground SEDs and provided the requirements on absolute $\alpha_f$ calibration are met. However, while those results are promising and seem robust against complex foregrounds and systematic effects, they also show the limits of the simple spatially constant foreground SED used. Indeed, results from Phase 2 onward show clear biases in $\alpha_f$ at low frequencies caused by the foregrounds' complexity. The estimation on $\beta$ is not significantly biased thanks to the compensation by precise $\alpha_f$ calibration at higher frequencies. Assessing the stability of systematic, foreground, and cosmological parameters across different sky fractions can serve as a robustness test against spatial variations in future uses of this pipeline. However, this demonstrates the need for an extension to the method which includes spatially variable foreground parameters to handle more complex foregrounds in the future. In particular, since J23 is based on *FGBuster*, a multi-patch like extension [106] constitute a natural next step. Indeed, it has been used in ref. [14] for *LiteBIRD* forecasts as well as for the estimation of the component separated maps used for the *D*-estimator and the pixel-based estimator in this paper. Other extensions capable of handling more complex foregrounds and systematic effects can be added to parametric methods in general although at the price of higher statistical uncertainty.

# 5 Discussion about the impact of foreground contamination

Our current understanding of Galactic emission is inherently limited by the data presently available, and a comprehensive grasp of the sky's true complexity will remain elusive until new data emerges. See ref. [108] and references therein, for example refs. [69, 109, 110], where new models of polarised Galactic dust and synchrotron emission at CMB frequencies
have been developed drawing on the latest observational constraints. We acknowledge the limitations of our simulations, particularly concerning the foreground models included since *pysm's* d0 and d1 models do not account for the line-of-sight variations in the dust SED, which can induce frequency-dependent changes in the polarisation angle as has been observed in *Planck* data [66, 67]. This frequency decorrelation may impact component separation methods and the characterization of primordial signals, including the *EB* spectrum [68, 73]. Regarding the synchrotron models s0 and s1, they show less spatial variability than what has been observed in ground-based measurements [111, 112]. Nevertheless, current constraints on the synchrotron *EB* signal are consistent with zero, suggesting a smaller impact compared to thermal dust [86, 87]. However, synchrotron emission remains relevant for accurately characterising polarisation angles. While our analysis does not explicitly model these complexities, our conclusions remain robust within the limitations of the foreground models considered. Future work should explore more sophisticated sky models (e.g., those of refs. [69–73]) as well as other instrumental effects due to, e.g., a non-ideal half-wave plate [113], to refine our understanding of these astrophysical and instrumental effects.

However, we expect that these additional foreground complexities will not substantially alter the estimated isotropic birefringence angle presented in this work. The main argument supporting this statement is the following: independently of their complexity, foreground models are expected to primarily impact the lowest multipoles. We have, however, shown for two out of the five pipelines presented here that our results are weakly dependent on the lowest $l$s. Regarding the *D*-estimator, we excluded the first 50 or 100 multipoles from the analysis in section 4.1.1 (see figure 6). Of course, this statement assumes that the component separation layer effectively handles foregrounds at multipoles larger than 100 and smaller than 750, thereby minimising their impact on our results within this range, which we believe to be reasonable. Similarly, although the MK estimator (see section 3.2.1), benefits from accessing the foreground signal at large angular scales to break the degeneracy between CB and polarisation angles, it can still tolerate the exclusion of the lowest multipoles if the foreground model is not reliable. We tested that an $l_{min} \approx 100$ cut leads to only <0.7 $\sigma$ shifts in $\beta$ compared to the results presented in section 4.2.1. Furthermore, we have seen that an incorrect modeling of the dust *EB* correlation could bias the CB angle measured with the MK estimator. However, ref. [70] showed that the difference between $C_{l}^{EB,dust}$ models can be partially absorbed by the extra degrees of freedom the $A$ amplitudes provide, resulting in <1 $\sigma$ shifts in $\beta$ as long as a reasonable template for *EB* is provided. Alternative ways to estimate the dust *EB* contribution from measurements of the magnetic misalignment angle also exist [53, 54].

In general, the component separation layer employed to provide the input for the methods in section 3.1 and all the techniques presented in section 3.2 can be extended to parametrize the foreground complexity. Of course, the increase in the number of parameters employed comes along with an increase in their statistical uncertainties which needs to be quantified. Ultimately, a robust assessment will require end-to-end simulations incorporating realistic noise properties, scanning strategy, and a representative level of instrumental systematics. Such simulations will be essential to draw firm conclusions regarding the impact of foreground complexity on the measurement of the birefringence angle which can be the subject for future analyses.
<table><caption>**Table 10.** Summary of the average performances of the five pipelines. Error bars show the standard deviation over simulations.</caption><thead><tr><th>Pipeline</th><th>Phase 1<br>β [×10<sup>-2</sup> deg]</th><th>Phase 2<br>β [×10<sup>-2</sup> deg]</th><th>Phase 3<br>β [×10<sup>-2</sup> deg]</th><th>Phase 4<br>β [×10<sup>-2</sup> deg]</th></tr></thead><tbody><tr><td>D-estimator</td><td>0.2 ± 1.0</td><td>0.1 ± 1.0</td><td>0.4 ± 4.6</td><td>29.7 ± 4.7</td></tr><tr><td>Pixel-based estimator</td><td>-0.3 ± 3.0</td><td>-0.5 ± 3.7</td><td>0.2 ± 4.1</td><td>28.1 ± 4.0</td></tr><tr><td>Template-based MK</td><td>-0.2 ± 1.5</td><td>-0.2 ± 1.7</td><td>-0.3 ± 2.1</td><td>29.3 ± 2.2</td></tr><tr><td>CAB-SeCRET</td><td>-0.4 ± 3.4</td><td>1.4 ± 3.4</td><td>1.7 ± 3.5</td><td>29.4 ± 3.6</td></tr><tr><td>J23</td><td>-0.2 ± 5.1</td><td>-1.0 ± 5.5</td><td>-1.8 ± 5.6</td><td>28.2 ± 5.6</td></tr></tbody></table>

<figure><img src="image_25.png" /><figcaption>**Figure 22.** Summary of this paper's results, showing the constraints derived from the five pipelines against the 4 phases of simulations.</figcaption></figure>

## 6 Conclusion

We generated four sets of 100 simulations, each containing CMB, foregrounds, and noise, with increasing complexity, as summarised in table 1. These simulations include cases with constant and spatially varying spectral indices of the foregrounds, with and without instrumental $\alpha_f$ rotation, and with and without a non-zero isotropic CB (see section 2 for details). Using five different pipelines, we recovered the mean $\beta$ summarised in table 10 and shown in figure 22.

Our results show that the constraints on $\beta$ are consistent across the four sets of simulations and all estimators, without significant bias from the instrumental or astrophysical systematic effects included in the simulations. More specifically, in most cases, the expected value for $\beta$ is recovered within 3 standard deviations of the mean. However, for certain specific pipelines and phases, we observe deviations ranging from 3 to almost 5 standard deviations of the mean, which correspond, in the worst case, to approximately 50% of the statistical uncertainty. Note that in Phases 3 and 4, the inclusion of non-zero $\alpha_f$ meets the requirements presented in ref. [60] and introduces both relative miscalibrations between frequency channels and a small global offset across all frequencies. Given the small value of the global offset, pipelines without control over $\alpha_f$ remain unbiased but reflect an increased uncertainty in $\beta$ due to the relative miscalibrations compared to Phases 1 and 2, which do not include these systematic effects.
The statistical efficiency with which we recover CB depends on the complexity of the simulations considered and the pipeline used. In general, we find that, for each pipeline, the error bar increases with the complexity of the simulations, with the uncertainty of $\beta$ ranging between 0.01° to 0.06° depending on the case considered. However, this effect is more pronounced for pipelines that ignore $\alpha_f$. Specifically, when we consider Phase 4, i.e., the case with $\beta = 0.3°$ motivated by recent observational hints [50], the five pipelines reach a detection with a significance from 5 to $13\sigma$. This is a remarkable result, since these estimates have been obtained considering realistic conditions containing the relevant systematic effects, see also discussion in section 5. To improve our grasp of these astrophysical and instrumental effects, it will be crucial for future research to investigate more advanced sky models (e.g., those by refs. [69–73]) and consider other instrumental influences, including those from a non-ideal half-wave plate [113].

Still focusing on Phase 4, we note that, among all the pipelines, the one providing the tightest constraints is the method based on the Minami-Komatsu approach (see sections 3.2.1 and 4.2.1). Compared to the D-estimator, which is also defined in harmonic space, the statistical power of the MK estimator comes from exploiting the cross-correlations between frequency bands instead of constraining $\beta$ from the cleaned CMB EB spectrum. It also achieves smaller error bars than the (β)CAB-SECRET and J23 pipelines as it does not marginalise over the uncertainties in foreground modeling. Although the MK, (β)CAB-SECRET, and J23 pipelines explicitly depend on our choice of foreground model (the D-estimator and stacking analysis indirectly do so through the previous component-separation step), they have so far adapted satisfactorily to the measurement of CB in realistic skies. We expect their performance to improve in the future as we continue to improve our knowledge of polarised foreground emission and implement clustering or foreground marginalisation techniques to adapt to more complex skies.

This MK technique, along with (β)CAB-SECRET and J23, demonstrates the capability to calibrate $\alpha_f$ and simultaneously detect CB under realistic conditions as those considered in the present analysis. The inclusion of information from instrumental calibration, whether from astrophysical sources [114] and laboratory or space-based instruments as proposed in refs. [115, 116], could improve the statistical efficiency of these methods and unambiguously break the degeneracy between isotropic CB and absolute miscalibration angles. It would also make the techniques that do not control for this systematic effect more robust.

Having five semi-independent pipelines aimed at estimating $\beta$ that produce consistent results is highly beneficial for ensuring the robustness of the analysis. This is because these pipelines respond differently to systematic effects, reducing the likelihood that a specific artefact affects all of them in the same way. Consequently, even though the pipelines have varying levels of statistical efficiency (see table 10), running all of them is essential for ensuring consistency and cross-validating the results. This approach significantly enhances the reliability of the analysis and the robustness of the results.

## Acknowledgments

This work is supported in Japan by ISAS/JAXA for Pre-Phase A2 studies, by the acceleration program of JAXA research and development directorate, by the World Premier Interna-
tional Research Center Initiative (WPI) of MEXT, by the JSPS Core-to-Core Program of A. Advanced Research Networks, and by JSPS KAKENHI Grant Numbers JP15H05891, JP17H01115, and JP17H01125. The Canadian contribution is supported by the Canadian Space Agency. The French *LiteBIRD* phase A contribution is supported by the Centre National d'Etudes Spatiale (CNES), by the Centre National de la Recherche Scientifique (CNRS), and by the Commissariat à l'Energie Atomique (CEA). The German participation in *LiteBIRD* is supported in part by the Excellence Cluster ORIGINS, which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy (Grant No. EXC-2094 - 390783311). The Italian *LiteBIRD* phase A contribution is supported by the Italian Space Agency (ASI Grants No. 2020-9-HH.0 and 2016-24-H.1-2018), the National Institute for Nuclear Physics (INFN) and the National Institute for Astrophysics (INAF). Norwegian participation in *LiteBIRD* is supported by the Research Council of Norway (Grant No. 263011) and has received funding from the European Research Council (ERC) under the Horizon 2020 Research and Innovation Programme (Grant agreement No. 772253 and 819478). The Spanish *LiteBIRD* phase A contribution is supported by MCIN/AEI/10.13039/501100011033, project refs. PID2019-110610RB-C21, PID2020-120514GB-I00, PID2022-139223OB-C21 (funded also by European Union NextGenerationEU/PRTR), and by MCIN/CDTI ICTP20210008 (funded also by EU FEDER funds). Funds that support contributions from Sweden come from the Swedish National Space Agency (SNSA/Rymdstyrelsen) and the Swedish Research Council (Reg. no. 2019-03959). The U.K. *LiteBIRD* contribution is supported by the U.K. Space Agency under grant reference ST/Y006003/1 - LiteBIRD U.K.: A major U.K. contribution to the LiteBIRD mission - Phase1 (March 25). The U.S. contribution is supported by NASA grant no. 80NSSC18K0132. We acknowledge funding from the SCIPOL project<fnref n="14" /> (PI: Josquin Errard, Grant agreement No. 101044073) and Marie Skłodowska-Curie program (Grant agreement No. 101007633) funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program. NB, AG, and PN acknowledge support by the MUR PRIN2022 Project BROWSEPOL: Beyond standaRd mOdel With coSmic microwavE background polarisation-2022EJNZ53 financed by the European Union-Next Generation EU. RS acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC). EdlH acknowledges support from IN2P3. This research used resources from the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility. We acknowledge the use of the **healpy** [117], **pysm** [62], **FGBuster** [74], **NaMaster** [92], **emcee** [94], **scipy** [97], ** numpy** [118], and **plotlib** [119] software packages. We thank Jiazheng Dou for carefully reviewing the equations in the appendix.

A Definition of the template-based Minami-Komatsu estimator

'In this work, we have extended the formalism presented in ref. [55] to include both synchrotron and dust templates. From eq. (3.10), we can build a Gaussian likelihood (introducing our eq. (3.10) into equation (A2) from ref. [55]) and simultaneously calculate the maximum-likelihood solution for all $\mathbf{x}_i = (\mathcal{A}^{\text{sync}}, \mathcal{A}^{\text{sync}\times\text{dust}}, \mathcal{A}^{\text{dust}}, \beta, \alpha_i)$ parameters in a semi-analytical

<footnote n="14">¹⁴https://scipol.in2p3.fr/</footnote>
iterative approach. The algorithm assumes that the $x_i$ parameters in the covariance matrix are known and fixed, starting at (1, 1, 1, 0, 0) values. We then obtain a system of linear equations to analytically calculate the maximum-likelihood solution by differentiating the log-likelihood with respect to each parameter. This first estimate is used to update the covariance matrix and recalculate a new best-fit solution, starting an iterative process that converges after only a few iterations.

For brevity, we introduce the $\sin(x) = s_x$ and $\cos(x) = c_x$ notation and the ‘o’, ‘s’, ‘d’, and ‘s × d’ shorthand for ‘obs’, ‘sync’, ‘dust’, and ‘sync×dust’ superscripts in the remainder of this section. Except for eq. (A.7), all foreground and CMB power spectra are implicitly convolved by the instrumental beam and pixel window functions associated with each frequency band. In the following equations, the ‘CMB’ angular power spectra correspond to the $\Lambda$CDM prediction for the CMB EE and BB spectra in the absence of CB.

The minimisation of the log-likelihood built from eq. (3.10) leads to a linear system of the form $\mathbf{A}_{mn}\mathbf{x}_m = \mathbf{b}_m$. The $\mathbf{A}_{mn}$ system matrix is symmetrical, with diagonal elements

$$ \mathbf{A}_{\mathrm{s,s}} = \sum_{\substack{b \\ i,j \neq i \\ p,q \neq p}} C_b^{E_i^s B_j^s} C_{ijpqb}^{-1} C_b^{E_p^s B_q^s}, $$

$$ \mathbf{A}_{\mathrm{s \times d, s \times d}} = \sum_{\substack{b \\ i,j \neq i \\ p,q \neq p}} \left[ C_b^{E_i^s B_j^d} C_{ijpqb}^{-1} C_b^{E_p^s B_q^d} + C_b^{E_i^d B_j^s} C_{ijpqb}^{-1} C_b^{E_p^d B_q^s} + 2C_b^{E_i^s B_j^d} C_{ijpqb}^{-1} C_b^{E_p^d B_q^s} \right], $$

$$ \mathbf{A}_{\mathrm{d,d}} = \sum_{\substack{b \\ i,j \neq i \\ p,q \neq p}} C_b^{E_i^d B_j^d} C_{ijpqb}^{-1} C_b^{E_p^d B_q^d}, $$

$$ \mathbf{A}_{\beta,\beta} = 4 \sum_{\substack{b \\ i,j \neq i \\ p,q \neq p}} \left( C_b^{E_i^{\mathrm{CMB}} E_j^{\mathrm{CMB}}} - C_b^{E_i^{\mathrm{CMB}} B_j^{\mathrm{CMB}}} \right) C_{ijpqb}^{-1} \left( C_b^{E_p^{\mathrm{CMB}} E_q^{\mathrm{CMB}}} - C_b^{E_p^{\mathrm{CMB}} B_q^{\mathrm{CMB}}} \right), \quad (A.1) $$

and off-diagonal terms

$$ \mathbf{A}_{\mathrm{s,s \times d}} = \sum_{\substack{b \\ i,j \neq i \\ p,q \neq p}} C_b^{E_i^s B_j^s} C_{ijpqb}^{-1} \left( C_b^{E_p^s B_q^d} + C_b^{E_p^d B_q^s} \right), $$

$$ \mathbf{A}_{\mathrm{s,d}} = \sum_{\substack{b \\ i,j \neq i \\ p,q \neq p}} C_b^{E_i^s B_j^s} C_{ijpqb}^{-1} C_b^{E_p^d B_q^d}, $$

$$ \mathbf{A}_{\mathrm{s},\beta} = 2 \sum_{\substack{b \\ i,j \neq i \\ p,q \neq p}} C_b^{E_i^s B_j^s} C_{ijpqb}^{-1} \left( C_b^{E_p^{\mathrm{CMB}} E_q^{\mathrm{CMB}}} - C_b^{E_p^{\mathrm{CMB}} B_q^{\mathrm{CMB}}} \right), $$

$$ \mathbf{A}_{\mathrm{s},\alpha_n} = 2 \sum_{\substack{b \\ i \neq n \\ p,q \neq p}} \left[ C_b^{E_i^o E_n^o} C_{inpqb}^{-1} C_b^{E_p^s B_q^s} - C_b^{B_n^o B_i^o} C_{nipqb}^{-1} C_b^{E_p^s B_q^s} \right], $$
$$ \mathbf{A}_{\text{s}\times\text{d},\text{d}} = \sum_{\substack{b \\ i,j \neq i \\ p,q \neq p}} C_b^{E_i^d B_j^d} C_{ijpqb}^{-1} (C_b^{E_p^s B_q^s} + C_b^{E_p^d B_q^s}), $$

$$ \mathbf{A}_{\text{s}\times\text{d},\beta} = 2 \sum_{\substack{b \\ i,j \neq i \\ p,q \neq p}} \left( C_b^{E_i^s B_j^d} + C_b^{E_i^d B_j^s} \right) C_{ijpqb}^{-1} \left( C_b^{E_p^{\text{CMB}} E_q^{\text{CMB}}} - C_b^{E_p^{\text{CMB}} B_q^{\text{CMB}}} \right), $$

$$ \mathbf{A}_{\text{s}\times\text{d},\alpha_n} = 2 \sum_{\substack{b \\ i \neq n \\ p,q \neq p}} \left( C_b^{E_i^o E_n^o} \mathbf{C}_{inpqb}^{-1} - C_b^{B_n^o B_i^o} \mathbf{C}_{nipqb}^{-1} \right) \left( C_b^{E_p^s B_q^d} + C_b^{E_p^s B_q^d} \right), $$

(A.2)

$$ \mathbf{A}_{\text{d},\beta} = 2 \sum_{\substack{b \\ i,j \neq i \\ p,q \neq p}} C_b^{E_i^d B_j^d} C_{ijpqb}^{-1} \left( C_b^{E_p^{\text{CMB}} E_q^{\text{CMB}}} - C_b^{B_p^{\text{CMB}} B_q^{\text{CMB}}} \right), $$

$$ \mathbf{A}_{\text{d},\alpha_n} = 2 \sum_{\substack{b \\ i \neq n \\ p,q \neq p}} \left[ C_b^{E_i^o E_n^o} \mathbf{C}_{inpqb}^{-1} C_b^{E_p^d B_q^d} - C_b^{B_n^o B_i^o} \mathbf{C}_{nipqb}^{-1} C_b^{E_p^d B_q^d} \right], $$

$$ \mathbf{A}_{\beta,\alpha_n} = 4 \sum_{\substack{b \\ i \neq n \\ p,q \neq p}} \left( C_b^{E_i^o E_n^o} \mathbf{C}_{inpqb}^{-1} - C_b^{B_n^o B_i^o} \mathbf{C}_{nipqb}^{-1} \right) \left( C_b^{E_p^{\text{CMB}} E_q^{\text{CMB}}} - C_b^{B_p^{\text{CMB}} B_q^{\text{CMB}}} \right) \quad (A.3) $$

$$ \mathbf{A}_{\alpha_m, \alpha_n} = 4 \sum_{\substack{b \\ i \neq n \\ j \neq m}} \left[ C_b^{E_i^o E_n^o} \mathbf{C}_{injmb}^{-1} C_b^{E_j^o E_m^o} + C_b^{B_n^o B_i^o} \mathbf{C}_{nimjb}^{-1} C_b^{B_m^o B_j^o} \right] \\ - 4 \sum_{\substack{b \\ i \neq n \\ j \neq m}} \left[ C_b^{B_n^o B_i^o} \mathbf{C}_{nijmb}^{-1} C_b^{E_j^o E_m^o} + C_b^{B_m^o B_j^o} \mathbf{C}_{mjinb}^{-1} C_b^{E_i^o E_n^o} \right]. $$

The elements of the independent term $\mathbf{b}_m$ are

$$ \mathbf{b}_s = \sum_{\substack{b \\ i,j \neq i \\ p,q \neq p}} C_b^{E_i^o B_j^o} C_{ijpqb}^{-1} C_b^{E_p^s B_q^s}, $$

$$ \mathbf{b}_{\text{s}\times\text{d}} = \sum_{\substack{b \\ i,j \neq i \\ p,q \neq p}} C_b^{E_i^o B_j^o} C_{ijpqb}^{-1} \left( C_b^{E_p^d B_q^s} + C_b^{E_p^s B_q^d} \right), $$

$$ \mathbf{b}_d = \sum_{\substack{b \\ i,j \neq i \\ p,q \neq p}} C_b^{E_i^o B_j^o} C_{ijpqb}^{-1} C_b^{E_p^d B_q^d}, \quad (A.4) $$

$$ \mathbf{b}_{\beta} = 2 \sum_{\substack{b \\ i,j \neq i \\ p,q \neq p}} C_b^{E_i^o B_j^o} C_{ijpqb}^{-1} \left( C_b^{E_p^{\text{CMB}} E_q^{\text{CMB}}} - C_b^{B_p^{\text{CMB}} B_q^{\text{CMB}}} \right), $$

$$ \mathbf{b}_{\alpha_m} = 2 \sum_{\substack{b \\ i \neq m \\ p,q \neq p}} \left[ C_b^{E_i^o E_m^o} \mathbf{C}_{impqb}^{-1} C_b^{E_p^o B_q^o} - C_b^{B_m^o B_i^o} \mathbf{C}_{mipqb}^{-1} C_b^{E_p^o B_q^o} \right]. $$
We exclude frequency auto-spectra ($i = j$ and $p = q$) from all summations on $\mathbf{A}_{mn}$ and $\mathbf{b}_m$ to avoid noise bias.

In our likelihood, we use an approximate covariance matrix derived from the observed signal on a realisation-by-realisation basis. For each combination of $i, j, p$, and $q$ frequency bands, the covariance can be divided into terms that depend only on the observed, foreground, and CMB spectra, and on their cross-correlations:

$$C_{ijpq\ell} = \frac{1}{(2\ell + 1)f_{\text{sky}}} \left[ C_{ijpq\ell}^{\text{fg}*\text{fg}} + C_{ijpq\ell}^{\text{fg}*\text{CMB}} - C_{ijpq\ell}^{\text{fg}*\text{fg}} + C_{ijpq\ell}^{\text{fg}*\text{CMB}} + C_{ijpq\ell}^{\text{fg}*\text{CMB}} \right], \quad (A.5)$$

where we have neglected any potential chance correlations between foreground and CMB signals ($C^{\text{fg}*\text{CMB}} = 0$) and included a sky-fraction factor ($f_{\text{sky}}$) to account for partial sky coverage. We make the following approximations when calculating the covariance matrix: we use observed power spectra instead of theoretical models; we neglect $\ell$-to-$\ell'$ correlations, even in the presence of a mask; we assume that the spherical harmonic coefficients of CMB, noise, and foreground signals are Gaussian and isotropic; and we discard the $C_{\ell}^{EB}$ terms suppressed by $\sin(x)^n$ with $n \ge 1$, since they rapidly fluctuate around zero and make the covariance matrix unstable.

The covariance of the observed signal is calculated as

$$\begin{aligned} C_{ijpq\ell}^{\text{fg}*\text{fg}} & \approx C_{\ell}^{E_i^o E_p^o} C_{\ell}^{B_j^o B_q^o} + C_{\ell}^{E_i^o B_q^o} C_{\ell}^{B_j^o E_p^o} \\ & \quad + \frac{s_{4\alpha_j} s_{4\alpha_q}}{[c_{4\alpha_i} + c_{4\alpha_j}][c_{4\alpha_p} + c_{4\alpha_q}]} \left( C_{\ell}^{E_i^o E_p^o} C_{\ell}^{E_j^o E_q^o} + C_{\ell}^{E_i^o E_q^o} C_{\ell}^{E_j^o E_p^o} \right) \quad (A.6) \\ & \quad + \frac{s_{4\alpha_i} s_{4\alpha_p}}{[c_{4\alpha_i} + c_{4\alpha_j}][c_{4\alpha_p} + c_{4\alpha_q}]} \left( C_{\ell}^{B_i^o B_p^o} C_{\ell}^{B_j^o B_q^o} + C_{\ell}^{B_i^o B_q^o} C_{\ell}^{B_j^o B_p^o} \right). \end{aligned}$$

The contribution of all CMB-related terms is

$$C_{ijpq\ell}^{\text{fg}*\text{CMB}} + C_{ijpq\ell}^{\text{fg}*\text{fg}} = - \frac{s_{4\beta}^2 b_{\ell}^i b_{\ell}^j b_{\ell}^p b_{\ell}^q \omega_{\text{pix},\ell}^A}{2c_{2\alpha_i+2\alpha_j}c_{2\alpha_p+2\alpha_q}} \left[ (C_{\ell}^{EE,\text{CMB}})^2 + (C_{\ell}^{BB,\text{CMB}})^2 \right], \quad (A.7)$$

where the combination of frequency bands is specified through the different beam and pixel window functions, $b_{\ell}^i$ and $\omega_{\text{pix},\ell}$, respectively.

In turn, the covariance associated with both foreground templates can be divided into four terms,

$$\begin{aligned} C_{ijpq\ell}^{\text{fg}*\text{fg}} = & \frac{4c_{2\alpha_i}c_{2\alpha_j}c_{2\alpha_p}c_{2\alpha_q}}{[c_{4\alpha_i} + c_{4\alpha_j}][c_{4\alpha_p} + c_{4\alpha_q}]} C_{ijpq\ell}^{\text{fg}_1} + \frac{4c_{2\alpha_i}c_{2\alpha_j}s_{2\alpha_p}s_{2\alpha_q}}{[c_{4\alpha_i} + c_{4\alpha_j}][c_{4\alpha_p} + c_{4\alpha_q}]} C_{ijpq\ell}^{\text{fg}_2} \\ & + \frac{4s_{2\alpha_i}s_{2\alpha_j}c_{2\alpha_p}c_{2\alpha_q}}{[c_{4\alpha_i} + c_{4\alpha_j}][c_{4\alpha_p} + c_{4\alpha_q}]} C_{ijpq\ell}^{\text{fg}_3} + \frac{4s_{2\alpha_i}s_{2\alpha_j}s_{2\alpha_p}s_{2\alpha_q}}{[c_{4\alpha_i} + c_{4\alpha_j}][c_{4\alpha_p} + c_{4\alpha_q}]} C_{ijpq\ell}^{\text{fg}_4}, \end{aligned} \quad (A.8)$$

which are calculated from the frequency cross-spectra of synchrotron and dust templates,
and their cross-correlations,

$$
\begin{equation}
\begin{aligned}
C_{ijpq\ell}^{\text{fg1}} \approx {}& (\mathcal{A}^s)^2 (C_{\ell}^{E_i^s E_p^s} C_{\ell}^{B_j^s B_q^s} + C_{\ell}^{E_i^s B_q^s} C_{\ell}^{B_j^s E_p^s}) + (\mathcal{A}^d)^2 (C_{\ell}^{E_i^d E_p^d} C_{\ell}^{B_j^d B_q^d} + C_{\ell}^{E_i^d B_q^d} C_{\ell}^{B_j^d E_p^d}) \\
& + (\mathcal{A}^{s \times d})^2 (C_{\ell}^{E_i^s E_p^s} C_{\ell}^{B_j^d B_q^d} + C_{\ell}^{E_i^s B_q^d} C_{\ell}^{B_j^d E_p^s} + C_{\ell}^{E_i^d E_p^d} C_{\ell}^{B_j^s B_q^s} + C_{\ell}^{E_i^d B_q^s} C_{\ell}^{B_j^s E_p^d}) \\
& + \mathcal{A}^s \mathcal{A}^d (C_{\ell}^{E_i^s E_p^d} C_{\ell}^{B_j^s B_q^d} + C_{\ell}^{E_i^s B_q^d} C_{\ell}^{B_j^s E_p^d} + C_{\ell}^{E_i^d E_p^s} C_{\ell}^{B_j^d B_q^s} + C_{\ell}^{E_i^d B_q^s} C_{\ell}^{B_j^d E_p^s}) \\
& + (\mathcal{A}^{s \times d})^2 (C_{\ell}^{E_i^s E_p^s} C_{\ell}^{B_j^d B_q^s} + C_{\ell}^{E_i^s B_q^s} C_{\ell}^{B_j^d E_p^s} + C_{\ell}^{E_i^d E_p^s} C_{\ell}^{B_j^s B_q^d} + C_{\ell}^{E_i^d B_q^d} C_{\ell}^{B_j^s E_p^s}) \\
& + \mathcal{A}^s \mathcal{A}^{s \times d} (C_{\ell}^{E_i^s E_p^d} C_{\ell}^{B_j^s B_q^s} + C_{\ell}^{E_i^s B_q^s} C_{\ell}^{B_j^d E_p^d} + C_{\ell}^{E_i^d E_p^s} C_{\ell}^{B_j^d B_q^s} + C_{\ell}^{E_i^d B_q^s} C_{\ell}^{B_j^d E_p^s}) \\
& + \mathcal{A}^s \mathcal{A}^{s \times d} (C_{\ell}^{E_i^d E_p^s} C_{\ell}^{B_j^d B_q^d} + C_{\ell}^{E_i^d B_q^d} C_{\ell}^{B_j^d E_p^s} + C_{\ell}^{E_i^s E_p^d} C_{\ell}^{B_j^d B_q^s} + C_{\ell}^{E_i^d B_q^s} C_{\ell}^{B_j^d E_p^s}) \\
& + \mathcal{A}^d \mathcal{A}^{s \times d} (C_{\ell}^{E_i^d E_p^d} C_{\ell}^{B_j^d B_q^s} + C_{\ell}^{E_i^d B_q^s} C_{\ell}^{B_j^d E_p^d} + C_{\ell}^{E_i^d E_p^d} C_{\ell}^{B_j^s B_q^d} + C_{\ell}^{E_i^d B_q^d} C_{\ell}^{B_j^s E_p^d}) \\
& + \mathcal{A}^d \mathcal{A}^{s \times d} (C_{\ell}^{E_i^d E_p^d} C_{\ell}^{B_j^d B_q^s} + C_{\ell}^{E_i^d B_q^s} C_{\ell}^{B_j^d E_p^d} + C_{\ell}^{E_i^d E_p^d} C_{\ell}^{B_j^s B_q^d} + C_{\ell}^{E_i^d B_q^d} C_{\ell}^{B_j^s E_p^d}),
\end{aligned}
\tag{A.9}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
C_{ijpq\ell}^{\text{fg2}} \approx {}& (\mathcal{A}^s)^2 C_{\ell}^{E_i^s E_q^s} C_{\ell}^{B_j^s B_p^s} + (\mathcal{A}^d)^2 C_{\ell}^{E_i^d E_q^d} C_{\ell}^{B_j^d B_p^d} + \mathcal{A}^s \mathcal{A}^d \left( C_{\ell}^{E_i^s E_q^d} C_{\ell}^{B_j^s B_p^d} + C_{\ell}^{E_i^d E_q^s} C_{\ell}^{B_j^d B_p^s} \right) \\
& + (\mathcal{A}^{s \times d})^2 \left( C_{\ell}^{E_i^s E_q^d} C_{\ell}^{B_j^d B_p^s} + C_{\ell}^{E_i^d E_q^s} C_{\ell}^{B_j^s B_p^d} + C_{\ell}^{E_i^s E_q^s} C_{\ell}^{B_j^d B_p^d} + C_{\ell}^{E_i^d E_q^d} C_{\ell}^{B_j^s B_p^s} \right) \\
& + \mathcal{A}^s \mathcal{A}^{s \times d} \left( C_{\ell}^{E_i^s E_q^d} C_{\ell}^{B_j^s B_p^s} + C_{\ell}^{E_i^d E_q^s} C_{\ell}^{B_j^d B_p^s} + C_{\ell}^{E_i^s E_q^s} C_{\ell}^{B_j^d B_p^d} + C_{\ell}^{E_i^d E_q^d} C_{\ell}^{B_j^s B_p^s} \right) \\
& + \mathcal{A}^d \mathcal{A}^{s \times d} \left( C_{\ell}^{E_i^d E_q^d} C_{\ell}^{B_j^d B_p^s} + C_{\ell}^{E_i^d B_q^s} C_{\ell}^{B_j^d B_p^d} + C_{\ell}^{E_i^d E_q^s} C_{\ell}^{B_j^d B_p^d} + C_{\ell}^{E_i^d E_q^d} C_{\ell}^{B_j^s B_p^d} \right),
\end{aligned}
\tag{A.10}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
C_{ijpq\ell}^{\text{fg3}} \approx {}& (\mathcal{A}^s)^2 C_{\ell}^{B_i^s B_q^s} C_{\ell}^{E_j^s E_p^s} + (\mathcal{A}^d)^2 C_{\ell}^{B_i^d B_q^d} C_{\ell}^{E_j^d E_p^d} + \mathcal{A}^s \mathcal{A}^d \left( C_{\ell}^{B_i^d B_q^s} C_{\ell}^{E_j^d E_p^s} + C_{\ell}^{E_j^s E_p^d} C_{\ell}^{B_i^s B_q^d} \right) \\
& + (\mathcal{A}^{s \times d})^2 \left( C_{\ell}^{B_i^s B_q^d} C_{\ell}^{E_j^d E_p^s} + C_{\ell}^{E_j^s E_p^d} C_{\ell}^{B_i^d B_q^s} + C_{\ell}^{B_i^d B_q^d} C_{\ell}^{E_j^s E_p^s} + C_{\ell}^{E_j^d E_p^d} C_{\ell}^{B_i^s B_q^s} \right) \\
& + \mathcal{A}^s \mathcal{A}^{s \times d} \left( C_{\ell}^{B_i^s B_q^s} C_{\ell}^{E_j^d E_p^s} + C_{\ell}^{E_j^s E_p^s} C_{\ell}^{B_i^s B_q^d} + C_{\ell}^{B_i^d B_q^s} C_{\ell}^{E_j^s E_p^s} + C_{\ell}^{E_j^s E_p^d} C_{\ell}^{B_i^s B_q^s} \right) \\
& + \mathcal{A}^d \mathcal{A}^{s \times d} \left( C_{\ell}^{B_i^s B_q^d} C_{\ell}^{E_j^d E_p^d} + C_{\ell}^{E_j^d E_p^s} C_{\ell}^{B_i^d B_q^d} + C_{\ell}^{B_i^d B_q^d} C_{\ell}^{E_j^s E_p^d} + C_{\ell}^{E_j^d E_p^d} C_{\ell}^{B_i^d B_q^s} \right),
\end{aligned}
\tag{A.11}
\end{equation}
$$

and

$$
\begin{equation}
\begin{aligned}
C_{ijpq\ell}^{\text{fg4}} \approx {}& (\mathcal{A}^s)^2 C_{\ell}^{B_i^s B_p^s} C_{\ell}^{E_j^s E_q^s} + (\mathcal{A}^d)^2 C_{\ell}^{B_i^d B_p^d} C_{\ell}^{E_j^d E_q^d} + \mathcal{A}^s \mathcal{A}^d \left( C_{\ell}^{B_i^s B_p^d} C_{\ell}^{E_j^s E_q^d} + C_{\ell}^{B_i^d B_p^s} C_{\ell}^{E_j^d E_q^s} \right) \\
& + (\mathcal{A}^{s \times d})^2 \left( C_{\ell}^{B_i^s B_p^s} C_{\ell}^{E_j^d E_q^d} + C_{\ell}^{B_i^d B_p^d} C_{\ell}^{E_j^s E_q^s} + C_{\ell}^{B_i^s B_p^d} C_{\ell}^{E_j^d E_q^s} + C_{\ell}^{B_i^d B_p^s} C_{\ell}^{E_j^s E_q^d} \right) \\
& + \mathcal{A}^s \mathcal{A}^{s \times d} \left( C_{\ell}^{B_i^s B_p^s} C_{\ell}^{E_j^s E_q^d} + C_{\ell}^{B_i^s B_p^s} C_{\ell}^{E_j^d E_q^s} + C_{\ell}^{B_i^s B_p^d} C_{\ell}^{E_j^s E_q^s} + C_{\ell}^{B_i^d B_p^s} C_{\ell}^{E_j^s E_q^s} \right) \\
& + \mathcal{A}^d \mathcal{A}^{s \times d} \left( C_{\ell}^{B_i^d B_p^s} C_{\ell}^{E_j^d E_q^d} + C_{\ell}^{B_i^s B_p^d} C_{\ell}^{E_j^d E_q^d} + C_{\ell}^{B_i^d B_p^d} C_{\ell}^{E_j^d E_q^s} + C_{\ell}^{B_i^d B_p^d} C_{\ell}^{E_j^s E_q^d} \right).
\end{aligned}
\tag{A.12}
\end{equation}
$$
Finally, the cross-correlation between the observations and both foreground templates is given as

$$ C_{ijpq\ell}^{\text{fg}*o} \approx \frac{2c_{2\alpha_i} c_{2\alpha_j}}{c_{4\alpha_i} + c_{4\alpha_j}} \left[ \begin{aligned} & A^s \left( C_{\ell}^{E_i^s E_p^o} C_{\ell}^{B_j^s B_q^o} + C_{\ell}^{E_i^s B_q^o} C_{\ell}^{B_j^s E_p^o} \right) + A^d \left( C_{\ell}^{E_i^d E_p^o} C_{\ell}^{B_j^d B_q^o} + C_{\ell}^{E_i^d B_q^o} C_{\ell}^{B_j^d E_p^o} \right) \end{aligned} \right] \\ & + \frac{2c_{2\alpha_i} c_{2\alpha_j}}{c_{4\alpha_i} + c_{4\alpha_j}} A^{s \times d} \left( C_{\ell}^{E_i^s E_p^o} C_{\ell}^{B_j^d B_q^o} + C_{\ell}^{E_i^s B_q^o} C_{\ell}^{B_j^d E_p^o} + C_{\ell}^{E_i^d E_p^o} C_{\ell}^{B_j^s B_q^o} + C_{\ell}^{E_i^d B_q^o} C_{\ell}^{B_j^s E_p^o} \right) \\ & + \frac{2c_{2\alpha_p} c_{2\alpha_q}}{c_{4\alpha_p} + c_{4\alpha_q}} \left[ A^s \left( C_{\ell}^{E_i^o E_p^o} C_{\ell}^{B_j^o B_q^s} + C_{\ell}^{E_i^o B_q^s} C_{\ell}^{B_j^o E_p^s} \right) + A^d \left( C_{\ell}^{E_i^o E_p^d} C_{\ell}^{B_j^o B_q^d} + C_{\ell}^{E_i^o B_q^d} C_{\ell}^{B_j^o E_p^d} \right) \right] \\ & + \frac{2c_{2\alpha_p} c_{2\alpha_q}}{c_{4\alpha_p} + c_{4\alpha_q}} A^{s \times d} \left( C_{\ell}^{E_i^o E_p^s} C_{\ell}^{B_j^o B_q^d} + C_{\ell}^{E_i^o B_q^d} C_{\ell}^{B_j^o E_p^s} + C_{\ell}^{E_i^o E_p^d} C_{\ell}^{B_j^o B_q^s} + C_{\ell}^{E_i^o B_q^s} C_{\ell}^{B_j^o E_p^d} \right) \\ & + \frac{2s_{2\alpha_i} s_{2\alpha_j}}{c_{4\alpha_i} + c_{4\alpha_j}} \left[ A^s C_{\ell}^{B_i^s B_q^o} C_{\ell}^{E_j^s E_p^o} + A^d C_{\ell}^{B_i^d B_q^o} C_{\ell}^{E_j^d E_p^o} + A^{s \times d} \left( C_{\ell}^{B_i^s B_q^o} C_{\ell}^{E_j^d E_p^o} + C_{\ell}^{B_i^d B_q^o} C_{\ell}^{E_j^s E_p^o} \right) \right] \\ & + \frac{2s_{2\alpha_p} s_{2\alpha_q}}{c_{4\alpha_p} + c_{4\alpha_q}} \left[ A^s C_{\ell}^{E_i^o E_q^s} C_{\ell}^{B_j^o B_p^s} + A^d C_{\ell}^{E_i^o E_q^d} C_{\ell}^{B_j^o B_p^d} + A^{s \times d} \left( C_{\ell}^{E_i^o E_q^d} C_{\ell}^{B_j^o B_p^s} + C_{\ell}^{E_i^o B_q^s} C_{\ell}^{B_j^o B_p^d} \right) \right] \end{aligned} \right] \quad (A.13) $$

## References

[1] S.M. Carroll, G.B. Field and R. Jackiw, Limits on a Lorentz and Parity Violating Modification of Electrodynamics, Phys. Rev. D 41 (1990) 1231 [INSPIRE].

[2] S.M. Carroll and G.B. Field, The Einstein equivalence principle and the polarization of radio galaxies, Phys. Rev. D 43 (1991) 3789 [INSPIRE].

[3] D. Harari and P. Sikivie, Effects of a Nambu-Goldstone boson on the polarization of radio galaxies and the cosmic microwave background, Phys. Lett. B 289 (1992) 67 [INSPIRE].

[4] S.M. Carroll, Quintessence and the rest of the world, Phys. Rev. Lett. 81 (1998) 3067 [astro-ph/9806099] [INSPIRE].

[5] V.A. Kostelecky and M. Mewes, Lorentz-violating electrodynamics and the cosmic microwave background, Phys. Rev. Lett. 99 (2007) 011601 [astro-ph/0702379] [INSPIRE].

[6] K.R.S. Balaji, R.H. Brandenberger and D.A. Easson, Spectral dependence of CMB polarization and parity, JCAP 12 (2003) 008 [hep-ph/0310368] [INSPIRE].

[7] D.J.E. Marsh, Axion Cosmology, Phys. Rep. 643 (2016) 1 [arXiv:1510.07633] [INSPIRE].

[8] E.G.M. Ferreira, Ultra-light dark matter, Astron. Astrophys. Rev. 29 (2021) 7 [arXiv:2005.03254] [INSPIRE].

[9] M. Li and X. Zhang, Cosmological CPT violating effect on CMB polarization, Phys. Rev. D 78 (2008) 103516 [arXiv:0810.0403] [INSPIRE].

[10] M. Pospelov, A. Ritz, C. Skordis, A. Ritz and C. Skordis, Pseudoscalar perturbations and polarization of the cosmic microwave background, Phys. Rev. Lett. 103 (2009) 051302 [arXiv:0808.0673] [INSPIRE].

[11] N.W. Boggess et al., The COBE mission — Its design and performance two years after launch, Astrophys. J. 397 (1992) 420 [INSPIRE].

[12] WMAP collaboration, The Microwave Anisotropy Probe (MAP) mission, Astrophys. J. 583 (2003) 1 [astro-ph/0301158] [INSPIRE].
[13] PLANCK collaboration, Planck 2018 results. I. Overview and the cosmological legacy of Planck, Astron. Astrophys. 641 (2020) A1 [arXiv:1807.06205] [INSPIRE].

[14] LITEBIRD collaboration, Probing Cosmic Inflation with the LiteBIRD Cosmic Microwave Background Polarization Survey, Prog. Theor. Exp. Phys. 2023 (2023) 042F01 [arXiv:2202.02773] [INSPIRE].

[15] BICEP/KECK collaboration, BICEP/Keck XII: Constraints on axionlike polarization oscillations in the cosmic microwave background, Phys. Rev. D 103 (2021) 042002 [arXiv:2011.03483] [INSPIRE].

[16] BICEP/KECK collaboration, BICEP/Keck XIV: Improved constraints on axionlike polarization oscillations in the cosmic microwave background, Phys. Rev. D 105 (2022) 022006 [arXiv:2108.03316] [INSPIRE].

[17] SPT-3G collaboration, Searching for axionlike time-dependent cosmic birefringence with data from SPT-3G, Phys. Rev. D 106 (2022) 042011 [arXiv:2203.16567] [INSPIRE].

[18] POLARBEAR collaboration, Constraints on axionlike polarization oscillations in the cosmic microwave background with POLARBEAR, Phys. Rev. D 108 (2023) 043017 [arXiv:2303.08410] [INSPIRE].

[19] F. Finelli and M. Galaverni, Rotation of Linear Polarization Plane and Circular Polarization from Cosmological Pseudo-Scalar Fields, Phys. Rev. D 79 (2009) 063002 [arXiv:0802.4210] [INSPIRE].

[20] M.A. Fedderke, P.W. Graham and S. Rajendran, Axion Dark Matter Detection with CMB Polarization, Phys. Rev. D 100 (2019) 015040 [arXiv:1903.02666] [INSPIRE].

[21] M. Kamionkowski, How to De-Rotate the Cosmic Microwave Background Polarization, Phys. Rev. Lett. 102 (2009) 111302 [arXiv:0810.1286] [INSPIRE].

[22] D. Contreras, P. Boubel and D. Scott, Constraints on direction-dependent cosmic birefringence from Planck polarization data, JCAP 12 (2017) 046 [arXiv:1705.06387] [INSPIRE].

[23] SPT collaboration, Searching for Anisotropic Cosmic Birefringence with Polarization Data from SPTpol, Phys. Rev. D 102 (2020) 083504 [arXiv:2006.08061] [INSPIRE].

[24] T. Namikawa et al., Atacama Cosmology Telescope: Constraints on cosmic birefringence, Phys. Rev. D 101 (2020) 083527 [arXiv:2001.10465] [INSPIRE].

[25] A. Gruppuso, D. Molinari, P. Natoli and L. Pagano, Planck 2018 constraints on anisotropic birefringence and its cross-correlation with CMB anisotropy, JCAP 11 (2020) 066 [arXiv:2008.10334] [INSPIRE].

[26] M. Bortolami, M. Billi, A. Gruppuso, P. Natoli and L. Pagano, Planck constraints on cross-correlations between anisotropic cosmic birefringence and CMB polarization, JCAP 09 (2022) 075 [arXiv:2206.01635] [INSPIRE].

[27] G. Zagatti, M. Bortolami, A. Gruppuso, P. Natoli, L. Pagano and G. Fabbian, Planck constraints on cosmic birefringence and its cross-correlation with the CMB, JCAP 05 (2024) 034 [arXiv:2401.11973] [INSPIRE].

[28] P.A.R. Ade et al., BICEP/Keck. XVII. Line-of-sight Distortion Analysis: Estimates of Gravitational Lensing, Anisotropic Cosmic Birefringence, Patchy Reionization, and Systematic Errors, Astrophys. J. 949 (2023) 43 [arXiv:2210.08038] [INSPIRE].

[29] A. Lue, L.-M. Wang and M. Kamionkowski, Cosmological signature of new parity violating interactions, Phys. Rev. Lett. 83 (1999) 1506 [astro-ph/9812088] [INSPIRE].
[30] B. Feng, H. Li, M.-z. Li and X.-m. Zhang, *Gravitational leptogenesis and its signatures in CMB*, Phys. Lett. B 620 (2005) 27 [hep-ph/0406269] [INSPIRE].

[31] G.-C. Liu, S. Lee and K.-W. Ng, *Effect on cosmic microwave background polarization of coupling of quintessence to pseudoscalar formed from the electromagnetic field and its dual*, Phys. Rev. Lett. 97 (2006) 161303 [astro-ph/0606248] [INSPIRE].

[32] E. Komatsu, *New physics from the polarized light of the cosmic microwave background*, Nat. Rev. Phys. 4 (2022) 452 [arXiv:2202.13919] [INSPIRE].

[33] Y.D. Takahashi et al., *Characterization of the BICEP Telescope for High-Precision Cosmic Microwave Background Polarimetry*, Astrophys. J. 711 (2010) 1141 [arXiv:0906.4069] [INSPIRE].

[34] B.J. Koopman, *Detector Development and Polarization Analyses for the Atacama Cosmology Telescope*, Ph.D. Thesis, Cornell University, Ithaca, NY, U.S.A. (2018) [DOI:10.7298/947m-mf06] [INSPIRE].

[35] PLANCK collaboration, *Planck intermediate results. XLIX. Parity-violation constraints from polarization data*, Astron. Astrophys. 596 (2016) A110 [arXiv:1605.08633] [INSPIRE].

[36] B.D. Sherwin and T. Namikawa, *Cosmic birefringence tomography and calibration independence with reionization signals in the CMB*, Mon. Not. Roy. Astron. Soc. 520 (2023) 3298 [arXiv:2108.09287] [INSPIRE].

[37] H. Nakatsuka, T. Namikawa and E. Komatsu, *Is cosmic birefringence due to dark energy or dark matter? A tomographic approach*, Phys. Rev. D 105 (2022) 123509 [arXiv:2203.08560] [INSPIRE].

[38] M. Galaverni, F. Finelli and D. Paoletti, *Redshift evolution of cosmic birefringence in CMB anisotropies*, Phys. Rev. D 107 (2023) 083529 [arXiv:2301.07971] [INSPIRE].

[39] A. Greco, N. Bartolo and A. Gruppuso, *A new solution for the observed isotropic cosmic birefringence angle and its implications for the anisotropic counterpart through a Boltzmann approach*, JCAP 10 (2024) 028 [arXiv:2401.07079] [INSPIRE].

[40] J.R. Eskilt, L. Herold, E. Komatsu, K. Murai, T. Namikawa and F. Naokawa, *Constraints on Early Dark Energy from Isotropic Cosmic Birefringence*, Phys. Rev. Lett. 131 (2023) 121001 [arXiv:2303.15369] [INSPIRE].

[41] T. Fujita, K. Murai, H. Nakatsuka and S. Tsujikawa, *Detection of isotropic cosmic birefringence and its implications for axionlike particles including dark energy*, Phys. Rev. D 103 (2021) 043509 [arXiv:2011.11894] [INSPIRE].

[42] A. Greco, N. Bartolo and A. Gruppuso, *Cosmic birefringence: cross-spectra and cross-bispectra with CMB anisotropies*, JCAP 03 (2022) 050 [arXiv:2202.04584] [INSPIRE].

[43] A. Greco, N. Bartolo and A. Gruppuso, *Probing Axions through Tomography of Anisotropic Cosmic Birefringence*, JCAP 05 (2023) 026 [arXiv:2211.06380] [INSPIRE].

[44] L. Yin, J. Kochappan, T. Ghosh and B.-H. Lee, *Is cosmic birefringence model-dependent?*, JCAP 10 (2023) 007 [arXiv:2305.07937] [INSPIRE].

[45] J. Kochappan, L. Yin, B.-H. Lee and T. Ghosh, *Observational evidence for Early Dark Energy as a unified explanation for Cosmic Birefringence and the Hubble tension*, arXiv:2408.09521 [INSPIRE].
[46] Y. Minami, H. Ochi, K. Ichiki, N. Katayama, E. Komatsu and T. Matsumura, *Simultaneous determination of the cosmic birefringence and miscalibrated polarization angles from CMB experiments*, Prog. Theor. Exp. Phys. **2019** (2019) 083E02 [arXiv:1904.12440] [INSPIRE].

[47] Y. Minami, *Determination of miscalibrated polarization angles from observed cosmic microwave background and foreground EB power spectra: Application to partial-sky observation*, Prog. Theor. Exp. Phys. **2020** (2020) 063E01 [arXiv:2002.03572] [INSPIRE].

[48] Y. Minami and E. Komatsu, *Simultaneous determination of the cosmic birefringence and miscalibrated polarization angles II: Including cross frequency spectra*, Prog. Theor. Exp. Phys. **2020** (2020) 103E02 [arXiv:2006.15982] [INSPIRE].

[49] Y. Minami and E. Komatsu, *New Extraction of the Cosmic Birefringence from the Planck 2018 Polarization Data*, Phys. Rev. Lett. **125** (2020) 221301 [arXiv:2011.11254] [INSPIRE].

[50] P. Diego-Palazuelos et al., *Cosmic Birefringence from the Planck Data Release 4*, Phys. Rev. Lett. **128** (2022) 091302 [arXiv:2201.07682] [INSPIRE].

[51] J.R. Eskilt and E. Komatsu, *Improved constraints on cosmic birefringence from the WMAP and Planck cosmic microwave background polarization data*, Phys. Rev. D **106** (2022) 063503 [arXiv:2205.13962] [INSPIRE].

[52] J.R. Eskilt, *Frequency-dependent constraints on cosmic birefringence from the LFI and HFI Planck Data Release 4*, Astron. Astrophys. **662** (2022) A10 [arXiv:2201.13347] [INSPIRE].

[53] S.E. Clark, C.-G. Kim, J.C. Hill and B.S. Hensley, *The Origin of Parity Violation in Polarized Dust Emission and Implications for Cosmic Birefringence*, Astrophys. J. **919** (2021) 53 [arXiv:2105.00120] [INSPIRE].

[54] A.J. Cukierman, S.E. Clark and G. Halal, *Magnetic Misalignment of Interstellar Dust Filaments*, Astrophys. J. **946** (2023) 106 [arXiv:2208.07382] [INSPIRE].

[55] P. Diego-Palazuelos et al., *Robustness of cosmic birefringence measurement against Galactic foreground emission and instrumental systematics*, JCAP **01** (2023) 044 [arXiv:2210.07655] [INSPIRE].

[56] COSMOGLOBE collaboration, *COSMOGLOBE DR1 results — II. Constraints on isotropic cosmic birefringence from reprocessed WMAP and Planck LFI data*, Astron. Astrophys. **679** (2023) A144 [arXiv:2305.02268] [INSPIRE].

[57] T. Louis et al., *The Atacama Cosmology Telescope: DR6 Power Spectra, Likelihoods and ΛCDM Parameters*, arXiv:2503.14452 [INSPIRE].

[58] A. Gruppuso and S. di Serego Alighieri, *Cosmic Polarization Rotation from CMB data: A review for GR110*, Int. J. Mod. Phys. D **34** (2025) 2540008 [arXiv:2502.07743] [INSPIRE].

[59] K.M. Górski et al., *HEALPix — A Framework for high resolution discretization, and fast analysis of data distributed on the sphere*, Astrophys. J. **622** (2005) 759 [astro-ph/0409513] [INSPIRE].

[60] P. Vielva et al., *Polarization angle requirements for CMB B-mode experiments. Application to the LiteBIRD satellite*, JCAP **04** (2022) 029 [arXiv:2202.01324] [INSPIRE].

[61] PLANCK collaboration, *Planck 2018 results. VI. Cosmological parameters*, Astron. Astrophys. **641** (2020) A6 [Erratum ibid. **652** (2021) C4] [arXiv:1807.06209] [INSPIRE].

[62] A. Zonca, B. Thorne, N. Krachmalnicoff and J. Borrill, *The Python Sky Model 3 software*, J. Open Source Softw. **6** (2021) 3783 [arXiv:2108.01444] [INSPIRE].
[63] PLANCK collaboration, Planck 2018 results. IV. Diffuse component separation, Astron. Astrophys. 641 (2020) A4 [arXiv:1807.06208] [INSPIRE].

[64] LITEBIRD collaboration, In-flight polarization angle calibration for LiteBIRD: blind challenge and cosmological implications, JCAP 01 (2022) 039 [arXiv:2111.09140] [INSPIRE].

[65] K.M. Huffenberger, A. Rotti and D.C. Collins, The Power Spectra of Polarized, Dusty Filaments, Astrophys. J. 899 (2020) 31 [arXiv:1906.10052] [INSPIRE].

[66] V. Pelgrims et al., Evidence for Line-of-Sight Frequency Decorrelation of Polarized Dust Emission in Planck Data, Astron. Astrophys. 647 (2021) A16 [arXiv:2101.09291] [INSPIRE].

[67] A. Ritacco et al., Dust polarization spectral dependence from Planck HFI data — Turning point for cosmic microwave background polarization-foreground modeling, Astron. Astrophys. 670 (2023) A163 [arXiv:2206.07671] [INSPIRE].

[68] L. Vacher et al., Frequency dependence of the thermal dust E/B ratio and EB correlation: Insights from the spin-moment expansion, Astron. Astrophys. 672 (2023) A146 [arXiv:2210.14768] [INSPIRE].

[69] L. Vacher, J. Chluba, J. Aumont, A. Rotti and L. Montier, High precision modeling of polarized signals: Moment expansion method generalized to spin-2 fields, Astron. Astrophys. 669 (2023) A5 [arXiv:2205.01049] [INSPIRE].

[70] C. Hervías-Caimapo, A.J. Cukierman, P. Diego-Palazuelos, K.M. Huffenberger and S.E. Clark, Modeling parity-violating spectra in Galactic dust polarization with filaments and its applications to cosmic birefringence searches, Phys. Rev. D 111 (2025) 083532 [arXiv:2408.06214] [INSPIRE].

[71] F. Vansyngel et al., Statistical simulations of the dust foreground to cosmic microwave background polarization, Astron. Astrophys. 603 (2017) A62 [Erratum ibid. 618 (2018) C4] [arXiv:1611.02577] [INSPIRE].

[72] G. Martínez-Solaeche, A. Karakci and J. Delabrouille, A 3D model of polarized dust emission in the Milky Way, Mon. Not. Roy. Astron. Soc. 476 (2018) 1310 [arXiv:1706.04162] [INSPIRE].

[73] L. Vacher et al., How bad could it be? Modelling the 3D complexity of the polarised dust signal using moment expansion, Astron. Astrophys. 697 (2025) A212 [arXiv:2411.11649] [INSPIRE].

[74] R. Stompor, S.M. Leach, F. Stivoli and C. Baccigalupi, Maximum Likelihood algorithm for parametric component separation in CMB experiments, Mon. Not. Roy. Astron. Soc. 392 (2009) 216 [arXiv:0804.2645] [INSPIRE].

[75] R. Stompor, J. Errard and D. Poletti, Forecasting performance of CMB experiments in the presence of complex foreground contaminations, Phys. Rev. D 94 (2016) 083526 [arXiv:1609.03807] [INSPIRE].

[76] J. Errard, F. Stivoli and R. Stompor, Framework for performance forecasting and optimization of CMB B-mode observations in presence of astrophysical foregrounds, Phys. Rev. D 84 (2011) 063005 [arXiv:1105.3859] [INSPIRE].

[77] QUAD collaboration, Parity Violation Constraints Using Cosmic Microwave Background Polarization Spectra from 2006 and 2007 Observations by the QUAD Polarimeter, Phys. Rev. Lett. 102 (2009) 161302 [arXiv:0811.0618] [INSPIRE].

[78] A. Gruppuso, G. Maggio, D. Molinari and P. Natoli, A note on the birefringence angle estimation in CMB data analysis, JCAP 05 (2016) 020 [arXiv:1604.05202] [INSPIRE].
[79] E. Komatsu et al., *Seven-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Interpretation*, Astrophys. J. Suppl. 192 (2011) 18 [arXiv:1001.4538] [INSPIRE].

[80] R.M. Sullivan, A. Abghari, P. Diego-Palazuelos, L.T. Hergt and D. Scott, *Planck PR4 (NPIPE) map-space cosmic birefringence*, JCAP 06 (2025) 025 [arXiv:2502.07654] [INSPIRE].

[81] D.L. Jow, D. Contreras, D. Scott and E.F. Bunn, *Taller in the saddle: constraining CMB physics using saddle points*, JCAP 03 (2019) 031 [arXiv:1811.05629] [INSPIRE].

[82] M. Kamionkowski, A. Kosowsky and A. Stebbins, *Statistics of cosmic microwave background polarization*, Phys. Rev. D 55 (1997) 7368 [astro-ph/9611125] [INSPIRE].

[83] A.P.S. Yadav, M. Su and M. Zaldarriaga, *Primordial B-mode Diagnostics and Self Calibrating the CMB Polarization*, Phys. Rev. D 81 (2010) 063512 [arXiv:0912.3532] [INSPIRE].

[84] N.J. Miller, M. Shimon and B.G. Keating, *CMB Polarization Systematics Due to Beam Asymmetry: Impact on Cosmological Birefringence*, Phys. Rev. D 79 (2009) 103002 [arXiv:0903.1116] [INSPIRE].

[85] W. Hu, M.M. Hedman and M. Zaldarriaga, *Benchmark parameters for CMB polarization experiments*, Phys. Rev. D 67 (2003) 043004 [astro-ph/0210096] [INSPIRE].

[86] F.A. Martire, R.B. Barreiro and E. Martínez-González, *Characterization of the polarized synchrotron emission from Planck and WMAP data*, JCAP 04 (2022) 003 [arXiv:2110.12803] [INSPIRE].

[87] J.A. Rubiño-Martín et al., *QUIJOTE scientific results — IV. A northern sky survey in intensity and polarization at 10–20 GHz with the Multi-Frequency Instrument*, Mon. Not. Roy. Astron. Soc. 519 (2023) 3383 [arXiv:2301.05113] [INSPIRE].

[88] H.K. Eriksen, J.B. Jewell, C. Dickinson, A.J. Banday, K.M. Górski and C.R. Lawrence, *Joint Bayesian component separation and CMB power spectrum estimation*, Astrophys. J. 676 (2008) 10 [arXiv:0709.1058] [INSPIRE].

[89] E. de la Hoz, P. Vielva, R.B. Barreiro and E. Martínez-González, *On the detection of CMB B-modes from ground at low frequency*, JCAP 06 (2020) 006 [arXiv:2002.12206] [INSPIRE].

[90] B. Thorne, J. Dunkley, D. Alonso and S. Næss, *The Python Sky Model: software for simulating the Galactic microwave sky*, Mon. Not. Roy. Astron. Soc. 469 (2017) 2821 [arXiv:1608.02841] [INSPIRE].

[91] E. de la Hoz, P. Diego-Palazuelos, E. Martínez-González, P. Vielva, R.B. Barreiro and J.D. Bilbao-Ahedo, *Determination of polarization angles in CMB experiments and application to CMB component separation analyses*, JCAP 03 (2022) 032 [arXiv:2110.14328] [INSPIRE].

[92] D. Alonso, J. Sanchez and A. Slosar, *A unified pseudo-Cℓ framework*, Mon. Not. Roy. Astron. Soc. 484 (2019) 4127 [arXiv:1809.09603] [INSPIRE].

[93] J.D. Bilbao-Ahedo, R.B. Barreiro, P. Vielva, E. Martínez-González and D. Herranz, *ECLIPSE: a fast Quadratic Maximum Likelihood estimator for CMB intensity and polarization power spectra*, JCAP 07 (2021) 034 [arXiv:2104.08528] [INSPIRE].

[94] D. Foreman-Mackey, D.W. Hogg, D. Lang and J. Goodman, *emcee: The MCMC Hammer*, Publ. Astron. Soc. Pac. 125 (2013) 306 [arXiv:1202.3665] [INSPIRE].

[95] B. Jost, J. Errard and R. Stompor, *Characterizing cosmic birefringence in the presence of Galactic foregrounds and instrumental systematic effects*, Phys. Rev. D 108 (2023) 082005 [arXiv:2212.08007] [INSPIRE].
[96] A. Rizzieri, J. Errard and R. Stompor, *Validating a main beam treatment of parametric, pixel-based component separation in the context of CMB observations*, Phys. Rev. D 111 (2025) 043512 [arXiv:2409.12944] [INSPIRE].

[97] P. Virtanen et al., *SciPy 1.0 — Fundamental Algorithms for Scientific Computing in Python*, Nat. Methods 17 (2020) 261 [arXiv:1907.10121] [INSPIRE].

[98] M. Tegmark, A.N. Taylor and A.F. Heavens, *Karhunen-Loeve eigenvalue problems in cosmology: How should we tackle large data sets?*, Astrophys. J. 480 (1997) 22 [astro-ph/9603021] [INSPIRE].

[99] A. Abghari, R.M. Sullivan, L.T. Hergt and D. Scott, *Constraints on cosmic birefringence using E-mode polarisation*, JCAP 08 (2022) 033 [arXiv:2203.10733] [INSPIRE].

[100] PLANCK collaboration, *Planck 2015 results. X. Diffuse component separation: Foreground maps*, Astron. Astrophys. 594 (2016) A10 [arXiv:1502.01588] [INSPIRE].

[101] M. Remazeilles, C. Dickinson, A.J. Banday, M.A. Bigot-Sazy and T. Ghosh, *An improved source-subtracted and destriped 408 MHz all-sky map*, Mon. Not. Roy. Astron. Soc. 451 (2015) 4311 [arXiv:1411.3628] [INSPIRE].

[102] C.L. Bennett et al., *Nine-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Final Maps and Results*, Astrophys. J. Suppl. 208 (2013) 20 [arXiv:1212.5225] [INSPIRE].

[103] M.A. Miville-Deschênes et al., *Separation of anomalous and synchrotron emissions using WMAP polarization data*, Astron. Astrophys. 490 (2008) 1093 [arXiv:0802.3345] [INSPIRE].

[104] L. Verde et al., *First year Wilkinson Microwave Anisotropy Probe (WMAP) observations: Parameter estimation methodology*, Astrophys. J. Suppl. 148 (2003) 195 [astro-ph/0302218] [INSPIRE].

[105] B.T. Draine and B. Hensley, *Magnetic Nanoparticles in the Interstellar Medium: Emission Spectrum and Polarization*, Astrophys. J. 765 (2013) 159 [arXiv:1205.7021] [INSPIRE].

[106] J. Errard and R. Stompor, *Characterizing bias on large scale CMB B-modes after galactic foregrounds cleaning*, Phys. Rev. D 99 (2019) 043529 [arXiv:1811.00479] [INSPIRE].

[107] K. Wolz et al., *The Simons Observatory: pipeline comparison and validation for large-scale B-modes*, Astron. Astrophys. 686 (2024) A16 [arXiv:2302.04276] [INSPIRE].

[108] PAN-EXPERIMENT GALACTIC SCIENCE GROUP, *Full-sky Models of Galactic Microwave Emission and Polarization at Sub-arcminute Scales for the Python Sky Model*, arXiv:2502.20452 [INSPIRE].

[109] A.M. Meisner and D.P. Finkbeiner, *Two-component Thermal Dust Emission Model: Application to the Planck HFI Maps*, in *American Astronomical Society Meeting Abstracts* 224, American Astronomical Society (2014), p. 220.14.

[110] H. Liu, J.-R. Li and Y.-F. Cai, *Evaluation of the Single-component Thermal Dust Emission Model in Cosmic Microwave Background Experiments*, Astrophys. J. Suppl. 276 (2025) 45 [arXiv:2411.04543] [INSPIRE].

[111] N. Krachmalnicoff et al., *S-PASS view of polarized Galactic synchrotron at 2.3 GHz as a contaminant to CMB observations*, Astron. Astrophys. 618 (2018) A166 [arXiv:1802.01145] [INSPIRE].

[112] E. de la Hoz et al., *QUIJOTE scientific results — VIII. Diffuse polarized foregrounds from component separation with QUIJOTE-MFI*, Mon. Not. Roy. Astron. Soc. 519 (2023) 3504 [arXiv:2301.05117] [INSPIRE].
[113] M. Monelli et al., *Impact of half-wave plate systematics on the measurement of cosmic birefringence from CMB polarization*, JCAP 03 (2023) 034 [arXiv:2211.05685] [INSPIRE].

[114] J. Aumont, J.F. Macías-Pérez, A. Ritacco, N. Ponthieu and A. Mangilli, *Absolute calibration of the polarisation angle for future CMB B-mode experiments from current and future measurements of the Crab nebula*, Astron. Astrophys. 634 (2020) A100 [arXiv:1805.10475] [INSPIRE].

[115] F.J. Casas et al., *L2-CalSat: A Calibration Satellite for Ultra-Sensitive CMB Polarization Space Missions*, Sensors 21 (2021) 3361 [INSPIRE].

[116] A. Ritacco et al., *Absolute Reference for Microwave Polarization Experiments. The COSMOCal Project and its Proof of Concept*, Publ. Astron. Soc. Pac. 136 (2024) 115001 [arXiv:2405.12135] [INSPIRE].

[117] A. Zonca et al., *healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in Python*, J. Open Source Softw. 4 (2019) 1298 [INSPIRE].

[118] C.R. Harris et al., *Array programming with NumPy*, Nature 585 (2020) 357 [arXiv:2006.10256] [INSPIRE].

[119] J.D. Hunter, *Matplotlib: A 2D Graphics Environment*, Comput. Sci. Eng. 9 (2007) 90 [INSPIRE].